pln=function(x,alpha,pars) { C=1 a=0.035 alpha1=1-alpha alpha2=alpha mu=pars[1] sigma1=exp(pars[2]) sigma2=exp(pars[3]) snsq=mean((x-mean(x))^2) part1=sum(log(alpha1*dnorm(x,mu,sigma1)+alpha2*dnorm(x,mu,sigma2))) part2=C*log(1-abs(1-2*alpha2)) part3=-a*(snsq/sigma1^2+log(sigma1^2/snsq)) -a*(snsq/sigma2^2+log(sigma2^2/snsq)) tot=-(part1+part2+part3) if ( (is.na(tot)) | (tot > 1e+10) | is.nan(tot) ) { tot=1e+10 } return(tot) } est.initial=function(x,alpha) { out2=c() len=25 for (i in 1:len) { pars=c( runif(1,mean(x)-2,mean(x)+2),log(runif(2,sd(x)/10,10*sd(x))) ) f=optim(par=pars,fn=pln,x=x,alpha=alpha) out1=c(alpha,f\$par,-f\$value) out2=rbind(out2,out1) } index=(1:len)[(out2[,5]>=max(out2[,5]))][1] out=out2[index,] out } est.ite=function(x,par) { output=rep(0,3) C=1 a=0.035 alpha=par[1] alpha1=1-alpha alpha2=alpha mu=par[2] sigma1=exp(par[3]) sigma2=exp(par[4]) output[1]=-pln(x,alpha,c(mu,log(sigma1),log(sigma2))) snsq=mean((x-mean(x))^2) for(i in 2:3) { pdf=alpha1*dnorm(x,mu,sigma1)+alpha2*dnorm(x,mu,sigma2) w1=alpha1*dnorm(x,mu,sigma1)/pdf w2=alpha2*dnorm(x,mu,sigma2)/pdf mu=( sigma2^2*sum(w1*x)+sigma1^2*sum(w2*x) )/( sigma2^2*sum(w1)+sigma1^2*sum(w2) ) sigma1=sqrt( (sum(w1*(x-mu)^2)+2*a*snsq)/(sum(w1)+2*a) ) sigma2=sqrt( (sum(w2*(x-mu)^2)+2*a*snsq)/(sum(w2)+2*a) ) if ( sum(w2)/(sum(w1)+sum(w2)) <= 0.5 ) { alpha=min( (sum(w2)+C)/(sum(w1)+sum(w2)+C), 0.5 ) } else { alpha=max( sum(w2)/(sum(w1)+sum(w2)+C), 0.5 ) } alpha1=1-alpha alpha2=alpha output[i]=-pln(x,alpha,c(mu,log(sigma1),log(sigma2))) } output } EMScale=function(x) { outem=c() for ( alpha in c(0.1,0.3,0.5) ) { initial=est.initial(x,alpha) pln.alt=est.ite(x,initial) outem=rbind(outem,pln.alt) } mu0=mean(x) sigma0=sqrt( mean( (x-mu0)^2 ) ) pln.null=-pln(x,0.5,c(mu0,log(sigma0),log(sigma0))) emstat=2*(apply(outem,2,max)-pln.null) round(as.numeric(emstat),4) }