myL <- function(p0,p1,n0,n1,n2) (p0^n0)*(p1^n1)*((1-p0-p1)^n2);  
myZ <- function(n,a,hatp0,hatp1) sqrt(n)*(4*(a-1)*hatp0+(6*a-4)*hatp1+(3-4*a))/sqrt(3-8*a+6*(a^2));

myCLRT <- function(n0,n1,n2){

n = n0+n1+n2
hatp0 = n0/n
hatp1 = n1/n

if ( (hatp0 <= 0.25) & (2*hatp0 <= hatp1) & (hatp1 <= 0.5) ) { 
CLRT = 2*(log(myL(hatp0,hatp1,n0,n1,n2))-log(myL(.25,.5,n0,n1,n2)) )
}

if (hatp1 > 0.5) { 
CLRT = ifelse(n0/(2*(n0+n2)) > 0.25, 0, 2*(log(myL(n0/(2*(n0+n2)),.5,n0,n1,n2))-log(myL(.25,.5,n0,n1,n2))) )

}

if (2*hatp0 > hatp1) {
CLRT = ifelse( (n1+n2)/(3*n) > 0.25, 0, 2*(log(myL((n1+n2)/(3*n),2*(n1+n2)/(3*n),n0,n1,n2))-log(myL(.25,.5,n0,n1,n2))) )
}

return(as.numeric(CLRT))

}


linkage <- function(n0,n1,n2,Nbootstrap){

n = n0+n1+n2;

hatp0 = n0/n;
hatp1 = n1/n;

st1 = 2*log(myL(hatp0,hatp1,n0,n1,n2)/myL(.25,.5,n0,n1,n2));
st2 = sqrt(2*n)*(1-2*hatp0-hatp1);
st3 = 4*sqrt(n/3)*(.75-hatp0-hatp1);
st4 = myZ(n,0.355,hatp0,hatp1);
st5 = max(myZ(n,0,hatp0,hatp1),myZ(n,0.5,hatp0,hatp1));
st6 = max(myZ(n,0,hatp0,hatp1),myZ(n,0.355,hatp0,hatp1),myZ(n,0.5,hatp0,hatp1));
st7 = myCLRT(n0,n1,n2);

# p1 = 1-pchisq(st1,2);
p2 = 1-pnorm(st2);
p3 = 1-pnorm(st3);
p4 = 1-pnorm(st4);

p1 = 0;
p5 = 0;
p6 = 0;
p7 = 0;

for(i in 1:Nbootstrap){

data = rmultinom(1,n,prob = c(0.25,0.5,0.25));

n0 = data[1];
n1 = data[2];
n2 = data[3];

hatp0 = n0/n;
hatp1 = n1/n;

p1 = p1 + I(2*log(myL(hatp0,hatp1,n0,n1,n2)/myL(.25,.5,n0,n1,n2)) > st1);
p5 = p5 + I(max(myZ(n,0,hatp0,hatp1),myZ(n,0.5,hatp0,hatp1)) > st5 );
p6 = p6 + I(max(myZ(n,0,hatp0,hatp1),myZ(n,0.355,hatp0,hatp1),myZ(n,0.5,hatp0,hatp1)) > st6);
p7 = p7 + I(myCLRT(n0,n1,n2) > st7);

}

out1 = c(round(st1,4),round(p1/Nbootstrap,4));
out2 = c(round(st2,4),round(p2,4));
out3 = c(round(st3,4),round(p3,4));
out4 = c(round(st4,4),round(p4,4));
out5 = c(round(st5,4),round(p5/Nbootstrap,4));
out6 = c(round(st6,4),round(p6/Nbootstrap,4));
out7 = c(round(st7,4),round(p7/Nbootstrap,4));

return(list(LRT = as.numeric(out1), mean = as.numeric(out2), 
            prop = as.numeric(out3), CLRT = as.numeric(out7), MERT = as.numeric(out4),
            M2 = as.numeric(out5), M3 = as.numeric(out6))); 

}


linkage(n0=22,n1=45,n2=33,Nbootstrap=10000)