f_full =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]
M <- sum(peddata[16:22])
lambdas <- c(exp(u6),exp(u5),exp(u5+b1),exp(u5),exp(u5+b1),
	exp(u4),exp(u4+b1+log(2)),exp(u4+b2),exp(u3+b1),exp(u3+b1),
	exp(u2+b1),exp(u2+b2),exp(u2+b1),exp(u2+b2),exp(u1+b2))
lambdamiss <- c(exp(u1+b2)+exp(u2+b2),exp(u2+b1)+exp(u3+b1),
	exp(u2+b2)+exp(u4+b2),exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1),
	exp(u3+b1)+exp(u5+b1),exp(u4)+exp(u5),exp(u5)+exp(u6))
sum(log(dpois(peddata[1:15],lambdas)))+   
(M>0)*sum(log(dpois(peddata[16:22],lambdamiss)))
}

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

## Second order derivative ##
d2f_full = function(pars){
u1 <- pars[1]
u2 <- pars[2]	
u3 <- pars[3]
u4 <- pars[4]	
u5 <- pars[5]
u6 <- pars[6]	
b1 <- pars[7]
b2 <- pars[8]	
d2f = NULL 
d2f = matrix(c(
#d2l/du1dx
-2*exp(u1+b2)+M22*exp(u1+u2+2*b2)/(exp(u1+b2)+exp(u2+b2))^2,
-M22*exp(u1+b2)*exp(u2+b2)/(exp(u1+b2)+exp(u2+b2))^2,0,0,0,0,
0,-2*exp(u1+b2),
#d2l/du2dx
-M22*exp(u1+b2)*exp(u2+b2)/(exp(u1+b2)+exp(u2+b2))^2,
-4*exp(u2+b1)-4*exp(u2+b2)+
M11*exp(u2+u4+u5+3*b1+log(2))/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2+
M12*exp(u2+u4+2*b2)/(exp(u2+b2)+exp(u4+b2))^2+
M21*exp(u2+u3+2*b1)/(exp(u2+b1)+exp(u3+b1))^2+
M22*exp(u1+u2+2*b2)/(exp(u1+b2)+exp(u2+b2))^2,
-M21*exp(u2+b1)*exp(u3+b1)/(exp(u2+b1)+exp(u3+b1))^2,
-M11*exp(u2+b1)*exp(u4+b1+log(2))/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2-
M12*exp(u2+b2)*exp(u4+b2)/(exp(u2+b2)+exp(u4+b2))^2,
-M11*exp(u2+b1)*exp(u5+b1)/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2,0,
-4*exp(u2+b1),-4*exp(u2+b2),
#d2l/du3dx
0,-M21*exp(u2+b1)*exp(u3+b1)/(exp(u2+b1)+exp(u3+b1))^2,
-4*exp(u3+b1)+M21*exp(u2+b1)*exp(u3+b1)/(exp(u2+b1)+exp(u3+b1))^2+
M01*exp(u3+b1)*exp(u5+b1)/(exp(u3+b1)+exp(u5+b1))^2,
0,-M01*exp(u3+b1)*exp(u5+b1)/(exp(u3+b1)+exp(u5+b1))^2,0,
-4*exp(u3+b1),0,
#d2l/du4dx
0,-M11*exp(u2+b1)*exp(u4+b1+log(2))/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2-M12*exp(u2+b2)*exp(u4+b2)/(exp(u2+b2)+exp(u4+b2))^2,
0,-2*exp(u4)-2*exp(u4+b1+log(2))-2*exp(u4+b2)+
M11*(exp(u2+b1)+exp(u5+b1))*exp(u4+b1+log(2))/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2+
M12*exp(u2+b2)*exp(u4+b2)/(exp(u2+b2)+exp(u4+b2))^2+
M10*exp(u4+u5)/(exp(u4)+exp(u5))^2,
-M11*exp(u4+b1+log(2))*exp(u5+b1)/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2-M10*exp(u4+u5)/(exp(u4)+exp(u5))^2,0,
-2*exp(u4+b1+log(2)),-2*exp(u4+b2),
#d2l/du5dx
0,-M11*exp(u2+b1)*exp(u5+b1)/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2,
-M01*exp(u3+b1)*exp(u5+b1)/(exp(u3+b1)+exp(u5+b1))^2,
-M11*exp(u4+b1+log(2))*exp(u5+b1)/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2-M10*exp(u4+u5)/(exp(u4)+exp(u5))^2,
-4*exp(u5)-4*exp(u5+b1)+
M11*(exp(u2+b1)+exp(u4+b1+log(2)))*exp(u5+b1)/(exp(u2+b1)+exp(u4+b1+log(2))+exp(u5+b1))^2+
M00*exp(u5+u6)/(exp(u5)+exp(u6))^2+M10*exp(u4+u5)/(exp(u4)+exp(u5))^2+
M01*exp(u3+b1)*exp(u5+b1)/(exp(u3+b1)+exp(u5+b1))^2,
-M00*exp(u5+u6)/(exp(u5)+exp(u6))^2,
-4*exp(u5+b1),0,
#d2l/du6dx
0,0,0,0,-M00*exp(u5+u6)/(exp(u5)+exp(u6))^2,
-2*exp(u6)+M00*exp(u5+u6)/(exp(u5)+exp(u6))^2,0,0,
#d2l/dbeta1dx
0,-4*exp(u2+b1),-4*exp(u3+b1),-2*exp(u4+b1+log(2)),-4*exp(u5+b1),0,
-4*exp(u2+b1)-4*exp(u3+b1)-2*exp(u4+b1+log(2))-4*exp(u5+b1),0,
#d2l/dbeta2dx
-2*exp(u1+b2),-4*exp(u2+b2),0,-2*exp(u4+b2),0,0,
0,-2*exp(u4+b2)-4*exp(u2+b2)-2*exp(u1+b2)
),ncol=8,byrow=T)
d2f
}