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)
}