f_red_zip =function(pars,peddata){
u1 <- pars[1]
u2 <- pars[2]	
u3 <- pars[3]
u4 <- pars[4]	
u5 <- pars[5]
u6 <- pars[6]	
pa <- pars[7]
M <- sum(peddata[16:22])	
log((peddata[1]==0)*(1-(1-pa)^4)+(1-pa)^4*dpois(peddata[1],exp(u6)))+ 
log((peddata[2]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[2],exp(u5)))+ 
log((peddata[3]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[3],exp(u5)))+  
log((peddata[4]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[4],exp(u5)))+  
log((peddata[5]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[5],exp(u5)))+  
log((peddata[6]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[6],exp(u4)))+
log((peddata[7]==0)*(1-2*(pa^2)*(1-pa)^2)+(2*(pa^2)*(1-pa)^2)*dpois(peddata[7],exp(u4+log(2))))+
log((peddata[8]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[8],exp(u4)))+
log((peddata[9]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[9],exp(u3)))+
log((peddata[10]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[10],exp(u3)))+
log((peddata[11]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[11],exp(u2)))+
log((peddata[12]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[12],exp(u2)))+
log((peddata[13]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[13],exp(u2)))+
log((peddata[14]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[14],exp(u2)))+
log((peddata[15]==0)*(1-pa^4)+(pa^4)*dpois(peddata[15],exp(u1)))+ 
(M>0)*(log(dpois(peddata[16],exp(u1)+exp(u2))) +
log(dpois(peddata[17],exp(u2)+exp(u3))) +
log(dpois(peddata[18],exp(u2)+exp(u4))) +
log(dpois(peddata[19],exp(u2)+exp(u4+log(2))+exp(u5))) +
log(dpois(peddata[20],exp(u3)+exp(u5))) +
log(dpois(peddata[21],exp(u4)+exp(u5))) +
log(dpois(peddata[22],exp(u5)+exp(u6))))
}

## First order derivative ##
df_red_zip = function(pars,peddata){
u1 <- pars[1]
u2 <- pars[2]	
u3 <- pars[3]
u4 <- pars[4]	
u5 <- pars[5]
u6 <- pars[6]	
pa <- pars[7]
M <- sum(peddata[16:22])	
c(
# dl/du1
((pa^4)*(peddata[15]-exp(u1))*dpois(peddata[15],exp(u1)))/((peddata[15]==0)*(1-pa^4)+(pa^4)*dpois(peddata[15],exp(u1)))+(M>0)*(peddata[16]*exp(u1)/(exp(u1)+exp(u2))-exp(u1)),
# dl/du2
((pa^3)*(1-pa)*(peddata[11]-exp(u2))*dpois(peddata[11],exp(u2)))/((peddata[11]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[11],exp(u2)))+
((pa^3)*(1-pa)*(peddata[12]-exp(u2))*dpois(peddata[12],exp(u2)))/((peddata[12]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[12],exp(u2)))+
((pa^3)*(1-pa)*(peddata[13]-exp(u2))*dpois(peddata[13],exp(u2)))/((peddata[13]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[13],exp(u2)))+
((pa^3)*(1-pa)*(peddata[14]-exp(u2))*dpois(peddata[14],exp(u2)))/((peddata[14]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[14],exp(u2)))+
(M>0)*(peddata[19]*exp(u2)/(exp(u2)+exp(u4+log(2))+exp(u5))-exp(u2)+ 
peddata[18]*exp(u2)/(exp(u2)+exp(u4))-exp(u2)+
peddata[17]*exp(u2)/(exp(u2)+exp(u3))-exp(u2)+
peddata[16]*exp(u2)/(exp(u1)+exp(u2))-exp(u2)),
# dl/du3
(((pa^2)*(1-pa)^2)*(peddata[9]-exp(u3))*dpois(peddata[9],exp(u3)))/((peddata[9]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[9],exp(u3)))+
(((pa^2)*(1-pa)^2)*(peddata[10]-exp(u3))*dpois(peddata[10],exp(u3)))/((peddata[10]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[10],exp(u3)))+
(M>0)*(peddata[17]*exp(u3)/(exp(u2)+exp(u3))-exp(u3)+
peddata[20]*exp(u3)/(exp(u3)+exp(u5))-exp(u3)),
# dl/du4
(((pa^2)*(1-pa)^2)*(peddata[6]-exp(u4))*dpois(peddata[6],exp(u4)))/((peddata[6]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[6],exp(u4)))+
(((2*pa^2)*(1-pa)^2)*(peddata[7]-exp(u4+log(2)))*dpois(peddata[7],exp(u4+log(2))))/((peddata[7]==0)*(1-2*(pa^2)*(1-pa)^2)+(2*(pa^2)*(1-pa)^2)*dpois(peddata[7],exp(u4+log(2))))+
(((pa^2)*(1-pa)^2)*(peddata[8]-exp(u4))*dpois(peddata[8],exp(u4)))/((peddata[8]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[8],exp(u4)))+
(M>0)*(peddata[19]*exp(u4+log(2))/(exp(u2)+exp(u4+log(2))+exp(u5))-exp(u4+log(2))+
peddata[18]*exp(u4)/(exp(u2)+exp(u4))-exp(u4)+
peddata[21]*exp(u4)/(exp(u4)+exp(u5))-exp(u4)),
# dl/du5
((pa*(1-pa)^3)*(peddata[2]-exp(u5))*dpois(peddata[2],exp(u5)))/((peddata[2]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[2],exp(u5)))+
((pa*(1-pa)^3)*(peddata[3]-exp(u5))*dpois(peddata[3],exp(u5)))/((peddata[3]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[3],exp(u5)))+
((pa*(1-pa)^3)*(peddata[4]-exp(u5))*dpois(peddata[4],exp(u5)))/((peddata[4]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[4],exp(u5)))+
((pa*(1-pa)^3)*(peddata[5]-exp(u5))*dpois(peddata[5],exp(u5)))/((peddata[5]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[5],exp(u5)))+
(M>0)*(peddata[19]*exp(u5)/(exp(u2)+exp(u4+log(2))+exp(u5))-exp(u5)+
peddata[22]*exp(u5)/(exp(u5)+exp(u6))-exp(u5)+
peddata[21]*exp(u5)/(exp(u4)+exp(u5))-exp(u5)+
peddata[20]*exp(u5)/(exp(u3)+exp(u5))-exp(u5)),
# dl/du6
(((1-pa)^4)*(peddata[1]-exp(u6))*dpois(peddata[1],exp(u6)))/((peddata[1]==0)*(1-(1-pa)^4)+(1-pa)^4*dpois(peddata[1],exp(u6)))+(M>0)*(peddata[22]*exp(u6)/(exp(u5)+exp(u6))-exp(u6)),
# dl/dpa
4*((1-pa)^3)*((peddata[1]==0)-dpois(peddata[1],exp(u6)))/((peddata[1]==0)*(1-(1-pa)^4)+(1-pa)^4*dpois(peddata[1],exp(u6)))+
((1-pa)^2)*(4*pa-1)*((peddata[2]==0)-dpois(peddata[2],exp(u5)))/((peddata[2]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[2],exp(u5)))+
((1-pa)^2)*(4*pa-1)*((peddata[3]==0)-dpois(peddata[3],exp(u5)))/((peddata[3]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[3],exp(u5)))+
((1-pa)^2)*(4*pa-1)*((peddata[4]==0)-dpois(peddata[4],exp(u5)))/((peddata[4]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[4],exp(u5)))+
((1-pa)^2)*(4*pa-1)*((peddata[5]==0)-dpois(peddata[5],exp(u5)))/((peddata[5]==0)*(1-pa*(1-pa)^3)+pa*(1-pa)^3*dpois(peddata[5],exp(u5)))+
pa*(1-pa)*(3*pa-2)*((peddata[6]==0)-dpois(peddata[6],exp(u4)))/((peddata[6]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[6],exp(u4)))+
2*pa*(1-pa)*(3*pa-2)*((peddata[7]==0)-dpois(peddata[7],exp(u4+log(2))))/((peddata[7]==0)*(1-2*(pa^2)*(1-pa)^2)+(2*(pa^2)*(1-pa)^2)*dpois(peddata[7],exp(u4+log(2))))+
pa*(1-pa)*(3*pa-2)*((peddata[8]==0)-dpois(peddata[8],exp(u4)))/((peddata[8]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[8],exp(u4)))+
pa*(1-pa)*(3*pa-2)*((peddata[9]==0)-dpois(peddata[9],exp(u3)))/((peddata[9]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[9],exp(u3)))+
pa*(1-pa)*(3*pa-2)*((peddata[10]==0)-dpois(peddata[10],exp(u3)))/((peddata[10]==0)*(1-(pa^2)*(1-pa)^2)+((pa^2)*(1-pa)^2)*dpois(peddata[10],exp(u3)))+
(pa^2)*(4*pa-3)*((peddata[11]==0)-dpois(peddata[11],exp(u2)))/((peddata[11]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[11],exp(u2)))+
(pa^2)*(4*pa-3)*((peddata[12]==0)-dpois(peddata[12],exp(u2)))/((peddata[12]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[12],exp(u2)))+
(pa^2)*(4*pa-3)*((peddata[13]==0)-dpois(peddata[13],exp(u2)))/((peddata[13]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[13],exp(u2)))+
(pa^2)*(4*pa-3)*((peddata[14]==0)-dpois(peddata[14],exp(u2)))/((peddata[14]==0)*(1-(pa^3)*(1-pa))+(pa^3)*(1-pa)*dpois(peddata[14],exp(u2)))+
-4*(pa^3)*((peddata[15]==0)-dpois(peddata[15],exp(u1)))/((peddata[15]==0)*(1-pa^4)+(pa^4)*dpois(peddata[15],exp(u1)))
)
}