#Clarke's data
#n1 = 186, n2 = 38, n3 = 13 and n4 = 284

nvec<-c(186,38,13,284)

th<-rep(1/3,3)
alpha.vec<-numeric(3)

abstol<-1e-16
par.change<-Inf
istep<-0
while(par.change > abstol){

	thold<-th
	alpha.vec[1]<-2*nvec[1]*th[1]^2/(th[1]^2+2*th[1]*th[3])+
			  nvec[3]+2*nvec[1]*th[1]*th[3]/(th[1]^2+2*th[1]*th[3])
	alpha.vec[2]<-2*nvec[2]*th[2]^2/(th[2]^2+2*th[2]*th[3])+
			  nvec[3]+2*nvec[2]*th[2]*th[3]/(th[2]^2+2*th[2]*th[3])
	alpha.vec[3]<-2*nvec[1]*th[1]*th[3]/(th[1]^2+2*th[1]*th[3])+
			  2*nvec[4]+2*nvec[2]*th[2]*th[3]/(th[2]^2+2*th[2]*th[3])

	th<-alpha.vec/sum(alpha.vec)
	
	like.val<-sum(alpha.vec*log(th))	

	par.change<-sqrt(sum((th-thold)^2))
	istep<-istep+1
	print(c(istep,th,like.val))
}