theta=0.1
p1=1/2+theta/4
p2=1/4-theta/4
p3=1/4-theta/4
p4=theta/4
x<-rmultinom(1,size=100,prob=c(p1,p2,p3,p4))

thetainitial=0.60  ####initial value, that is \theta_{0}
it  <- 0  ####iterative count
del <- 1  ####iterative adjustment
thetaold=thetainitial  ###assign the initial value to theta_{1}
while(abs(del) > 0.000001 && (it <- it+1) < 20) ###Loop for 20 iterations and
                                                ###prespecified presicion
	{
	  u2=(x[1]*thetaold)/(2+thetaold)        ###Expectation step
	  thetanew=(u2+x[4])/(sum(x)-x[1]+u2)    ###Maximization step
	  del=thetanew-thetaold                  ###Calculate the difference between two iterations
	  thetaold=thetanew                      ###assign thetanew to thetaold for the recursions
	  cat(it, thetanew, "\n")                ###List the iterations
	}




# Generate data from an exponential with theta=2, and with the second
# experiment truncated at t=3. Note that R uses a form of the
# exponential in which the parameter is a multiplier; i.e., the R
# parameter is 1/theta. Set the seed, so computations are reproducible.

set.seed(4)
n<-100
m<-500
theta<-2
t<-3
x<-rexp(n,1/theta)
r<-min(which(sort(rexp(m,1/theta))>=3))-1

# We begin with theta=1.
# (Note that theta.k is set to theta.kp1 at the beginning of the loop.)
theta.k<-.01
theta.kp1<-1
# Do some preliminary computations
n.xbar<-sum(x)
# Then loop and test for convergence
it  <- 0  ####iterative count
del <- 1  ####iterative adjustment
while(abs(del) > 0.000001 && (it <- it+1) < 20) ###Loop for 20 iterations
                                                ###and prespecified presicion
	{
	  theta.kp1<-(n.xbar+(m-r)*(t+theta.k)+r*(theta.k-t*exp(-t/theta.k)/(1-exp(-t/theta.k))))/(n+m)
	  del<- theta.kp1-theta.k
	  theta.k<-theta.kp1
	  cat(it, theta.kp1, "\n")
	}




# Normal mixture. Generate data from normal mixture with w=0.7,
# mu_1=0, sigma^2_1=1, mu_2=1, sigma^2_2=2.
# Note that R uses sigma, rather than sigma^2 in rnorm.
# Set the seed, so computations are reproducible.

set.seed(4)
n<-300
w<-0.7
mu1<-0
sigma21<-1
mu2<-5
sigma22<-2
x<-ifelse(runif(n)<w, rnorm(n,mu1,sqrt(sigma21)),rnorm(n,mu2,sqrt(sigma22)))

# Initialize
theta.k<-.1
theta.kp1<-.5

it  <- 0  ####iterative count
del <- 1  ####iterative adjustment
while(abs(del) > 0.000001 && (it <- it+1) < 20) ###Loop for 20 iterations
                                                ###and prespecified presicion
	{
	  tmp<-theta.k*dnorm(x,mu1,sqrt(sigma21))
	  ehat.k<-tmp/(tmp+(1-theta.k)*dnorm(x,mu2,sqrt(sigma22)))
	  theta.kp1<-mean(ehat.k)
	  del<- theta.kp1-theta.k
	  theta.k<-theta.kp1
	  cat(it, theta.kp1, "\n")
	}