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

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