p_zip <- function(pars,MFCdata,M22,M21,M12,M11,M01,M10,M00){
b1 <- pars[1]
b2 <- pars[2]
pa <- pars[3]
M <- M22+M21+M12+M11+M01+M10+M00
p <- rep(0,15)
c <- 1 + 2*pa*(exp(b1)-1)+(pa^2)*(1-2*exp(b1)+exp(b2))
p[1] = pa^0 * (1-pa)^4 / c
p[2] = pa^1 * (1-pa)^3 / c
p[3] = pa^1 * (1-pa)^3 * exp(b1) / c
p[4] = pa^1 * (1-pa)^3 / c
p[5] = pa^1 * (1-pa)^3 * exp(b1) / c
p[6] = pa^2 * (1-pa)^2 / c
p[7] = pa^2 * (1-pa)^2 * exp(b1) * 2 / c
p[8] = pa^2 * (1-pa)^2 * exp(b2) / c
p[9] = pa^2 * (1-pa)^2 * exp(b1) / c
p[10] = pa^2 * (1-pa)^2 * exp(b1) / c
p[11] = pa^3 * (1-pa)^1 * exp(b1) / c
p[12] = pa^3 * (1-pa)^1 * exp(b2) / c
p[13] = pa^3 * (1-pa)^1 * exp(b1) / c
p[14] = pa^3 * (1-pa)^1 * exp(b2) / c
p[15] = pa^4 * (1-pa)^0 * exp(b2) / c
MFCdata[1]*log(p[1])+sum(MFCdata[c(2,4)])*log(p[2])+sum(MFCdata[c(3,5)])*log(p[3])+sum(MFCdata[6])*log(p[6])+sum(MFCdata[7])*log(p[7])+sum(MFCdata[8])*log(p[8])+sum(MFCdata[9:10])*log(p[9])+sum(MFCdata[c(11,13)])*log(p[11])+sum(MFCdata[c(12,14)])*log(p[12])+sum(MFCdata[15])*log(p[15])+(M>0)*(M00*log((1-pa)^3)+M10*log(pa*(1-pa)^2)+M01*log(exp(b1)*pa*(1-pa)^2)+M11*log(exp(b1)*pa*(1-pa))+M21*log(exp(b1)*(pa^2)*(1-pa))+M12*log(exp(b2)*(pa^2)*(1-pa))+M22*log(exp(b2)*(pa^3))-M*log(c))
	}
	
dp_zip <- function(pars,MFCdata,M22,M21,M12,M11,M01,M10,M00){
b1 <- pars[1]
b2 <- pars[2]
pa <- pars[3]
M <- M22+M21+M12+M11+M01+M10+M00
c <- 1 + 2*pa*(exp(b1)-1)+(pa^2)*(1-2*exp(b1)+exp(b2))
c(sum(MFCdata[c(3,5,7,9:11,13)])+(M>0)*(M01+M11+M21)-(sum(MFCdata[1:15])+M)*(2*pa*(1-pa)*exp(b1))/c,
sum(MFCdata[c(8,12,14,15)])+(M>0)*(M12+M22)-(sum(MFCdata[1:15])+M)*((pa^2)*exp(b2))/c,
(sum(MFCdata[2:5])+2*sum(MFCdata[6:10])+3*sum(MFCdata[11:14])+4*N222+(M>0)*(M01+M10+M11+2*(M12+M21)+3*M22))/(pa*(1-pa))-(4*sum(MFCdata[1:15])+3*(M00+M01+M10+M12+M21+M22)+2*M11)/(1-pa)-(2*(sum(MFCdata[1:15])+M)*(exp(b1)-1+pa*(1-2*exp(b1)+exp(b2))))/c)
	}