R Code Mixtures


y=c(rnorm(70),runif(130,-1.5,1.5))
hist(y,breaks=15)
N=length(y)

#make space
M=2000
z=matrix(nrow=M,ncol=N)
p=1:M

#start values
p[1]=.5
z[1,]=rbinom(N,1,p[1])

dny=dnorm(y)
duy=dunif(y,-1.5,1.5)

print(p)

for (m in 2:M)
{
p.unif=p[m-1]*duy
p.norm=(1-p[m-1])*dny
p.star=p.unif/(p.unif+p.norm)
z[m,]=rbinom(N,1,p.star)
p[m]=rbeta(1,1+sum(z[m,]),1+N-sum(z[m,]))
}
post.z=apply(z,2,mean)

####################################################

Research Methods 3020