Nothing
# Chapter 6 R commands
library(combinat)
data(datha)
datha=as.matrix(datha)
hist(datha,nclas=200,xlab="",xlim=c(min(datha),max(datha)),ylab="",prob=TRUE,main="")
S=readline(prompt="Type <Return> to continue : ")
# Section 6.3
n=length(datha)
meand=mean(datha)
vard=var(datha)
sdd=sqrt(vard)
dat=plotmix()$sample
S=readline(prompt="Type <Return> to continue : ")
omega=function(z,x,p)
{
n=length(x)
n1=sum(z==1)
n2=n-n1
if (n1==0) xbar1=0 else xbar1=sum((z==1)*x)/n1
if (n2==0) xbar2=0 else xbar2=sum((z==2)*x)/n2
ss1=sum((z==1)*(x-xbar1)^2)
ss2=sum((z==2)*(x-xbar2)^2)
sqrt((n1+.25)*(n2+.25))*p^n1*(1-p)^n2*
exp(-((n1+.25)*ss1+(n2+.25)*ss2)/2)*
exp(-(n1*xbar1^2+n2*xbar2)/8)
}
omega(z=sample(1:2,4,rep=TRUE),x=plotmix(n=4,plot=FALSE)$samp,p=.8)
S=readline(prompt="Type <Return> to continue : ")
# Section 6.4
dat=plotmix()$sample
simu=gibbsmean(0.7,dat)
points(simu,pch=".")
S=readline(prompt="Type <Return> to continue : ")
mix=list(k=3,mu=mean(datha),sig=var(datha))
simu=gibbsnorm(100,datha,mix)
hist(datha,prob=T,main="",xlab="",ylab="",nclass=100,col="wheat")
x=y=seq(min(datha),max(datha),length=150)
yy=matrix(0,ncol=150,nrow=100)
for (i in 1:150)
{
yy[,i]=apply(simu$p*dnorm(x[i],mean=simu$mu,sd=sqrt(simu$sig)),1,sum)
y[i]=mean(yy[,i])
}
for (t in 51:100) lines(x,yy[t,],col="gold")
lines(x,y,lwd=2.3,col="sienna2")
S=readline(prompt="Type <Return> to continue : ")
dat=plotmix()$sample
simean=hmmeantemp(dat,10000,var=1)
points(simean,pch=19,cex=.4)
S=readline(prompt="Type <Return> to continue : ")
# Section 6.5
# recycling former simu
par(mfrow=c(3,2),mar=c(4,4,1,1.5))
plot(simu$mu[,1],ylim=range(simu$mu),ylab=expression(mu[i]),xlab="n",type="l",col="gold4")
lines(simu$mu[,2],col="sienna3")
lines(simu$mu[,3],col="steelblue")
plot(simu$mu[,2],simu$p[,2],col="sienna3",
xlim=range(simu$mu),ylim=range(simu$p),
xlab=expression(mu[i]),ylab=expression(p[i]))
points(simu$mu[,3],simu$p[,3],col="steelblue")
points(simu$mu[,1],simu$p[,1],col="gold4")
plot(simu$p[,1],ylim=range(simu$p),ylab=expression(p[i]),xlab="n",type="l",col="gold4")
lines(simu$p[,2],col="sienna3")
lines(simu$p[,3],col="steelblue")
plot(simu$p[,2],simu$sig[,2],col="sienna3",
xlim=range(simu$p),ylim=range(simu$sig),
xlab=expression(p[i]),ylab=expression(sigma[i]^2))
points(simu$p[,3],simu$sig[,3],col="steelblue")
points(simu$p[,1],simu$sig[,1],col="gold4")
plot(simu$sig[,1],ylim=range(simu$sig),ylab=expression(sigma[i]^2),xlab="n",type="l",col="gold4")
lines(simu$sig[,2],col="sienna3")
lines(simu$sig[,3],col="steelblue")
plot(simu$sig[,2],simu$mu[,2],col="sienna3",
xlim=range(simu$sig),ylim=range(simu$mu),
xlab=expression(sigma[i]^2),ylab=expression(mu[i]))
points(simu$sig[,3],simu$mu[,3],col="steelblue")
points(simu$sig[,1],simu$mu[,1],col="gold4")
S=readline(prompt="Type <Return> to continue : ")
indimap=order(simu$lopost,decreasing=TRUE)[1]
map=list(mu=simu$mu[indimap,],
sig=simu$sig[indimap,],
p=simu$p[indimap,])
lili=alloc=matrix(0,length(datha),3)
for (t in 1:length(datha))
{
lili[t,]=map$p*dnorm(datha[t],mean=map$mu,
sd=sqrt(map$sig))
lili[t,]=lili[t,]/sum(lili[t,])
}
ormu=orsig=orp=matrix(0,ncol=3,nrow=1000)
perma=permn(3)
for (t in 1:1000)
{
entropies=rep(0,factorial(3))
for (j in 1:n)
{
alloc[j,]=simu$p[t,]*dnorm(datha[j],mean=simu$mu[t,],
sd=sqrt(simu$sig[t,]))
alloc[j,]=alloc[j,]/sum(alloc[j,])
for (i in 1:factorial(3))
entropies[i]=entropies[i]+
sum(lili[j,]*log(alloc[j,perma[[i]]]))
}
best=order(entropies,decreasing=TRUE)[1]
ormu[t,]=simu$mu[t,perma[[best]]]
orsig[t,]=simu$sig[t,perma[[best]]]
orp[t,]=simu$p[t,perma[[best]]]
}
# Section 6.7
mu1=2.5;mu2=0;p=.7;n=500;nl=50
pbar=1-p
u=runif(n)
sampl=rnorm(n)+(u<=p)*mu1+(u>p)*mu2
mu1=mu2=seq(min(sampl),max(sampl),.1)
mo1=mu1%*%t(rep(1,length(mu2)))
mo2=rep(1,length(mu2))%*%t(mu2)
ca1=-0.5*mo1*mo1
ca2=-0.5*mo2*mo2
like=0*mo1
for (i in 1:n)
like=like+log(p*exp(ca1+sampl[i]*mo1)+
pbar*exp(ca2+sampl[i]*mo2))
like=like+.1*(ca1+ca2)
par(mfrow=c(1,3),mar=c(4,4,1,1))
for (i in 1:3){
image(mu1,mu2,like,xlab=expression(mu[1]),
ylab=expression(mu[2]),col=heat.colors(250))
contour(mu1,mu2,like,levels=seq(min(like),max(like),length=nl),
add=TRUE,drawlabels=FALSE)
simu=hmmeantemp(dat,10000,var=0.1,alpha=1/(10)^(i-1))
points(simu,pch=19,cex=.4)
}
S=readline(prompt="Type <Return> to continue : ")
dat=plotmix()$sample
simu=hmmeantemp(dat,10000)
points(simu,pch=19,cex=.4,col="sienna")
simu=hmmeantemp(dat,10000,alpha=0.1)
points(simu,pch=19,cex=.4,col="steelblue4")
simu=hmmeantemp(dat,10000,alpha=0.01)
points(simu,pch=19,cex=.4)
S=readline(prompt="Type <Return> to continue : ")
# Section 6.8
########## k=2
k=2
mu=rnorm(k,mean=meand,sd=sdd/k)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(1,k))
mix=list(k=k,p=p,mu=mu,sig=sig)
resim2=gibbs(200,datha,mix)
lulu=order(resim2$lolik)[200]
lnum1=resim2$lolik[lulu]
lnum2=sum(dnorm(resim2$mu[lulu,],mean=meand,sd=resim2$sig[lulu,],log=TRUE)+
dgamma(1/resim2$sig[lulu,],10,var(datha),log=TRUE)-2*log(resim2$sig[lulu,]))+
sum((rep(0.5,2)-1)*log(resim2$p[lulu,]))+
lgamma(sum(rep(0.5,2)))-sum(lgamma(rep(0.5,2)))
lchibapprox2=lnum1+lnum2-log(resim2$deno)
lchibapprox2
S=readline(prompt="Type <Return> to continue : ")
########## k=3
k=3
mu=rnorm(k,mean=meand,sd=sdd/k)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(1,k))
mix=list(k=k,p=p,mu=mu,sig=sig)
resim3=gibbs(200,datha,mix)
lulu=order(resim3$lolik)[200]
lnum1=resim3$lolik[lulu]
lnum2=sum(dnorm(resim3$mu[lulu,],mean=meand,sd=resim3$sig[lulu,],log=TRUE)+
dgamma(resim3$sig[lulu,],10,vard,log=TRUE)-2*log(resim3$sig[lulu,]))+
sum((rep(0.5,mix$k)-1)*log(resim3$p[lulu,]))+
lgamma(sum(rep(0.5,mix$k)))-sum(lgamma(rep(0.5,mix$k)))
lchibapprox3=lnum1+lnum2-log(resim3$deno)
S=readline(prompt="Type <Return> to continue : ")
########## k=4
k=4
mu=rnorm(k,mean=meand,sd=sdd/k)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(1,k))
mix=list(k=k,p=p,mu=mu,sig=sig)
resim4=gibbs(200,datha,mix)
lulu=order(resim4$lolik)[200]
lnum1=resim4$lolik[lulu]
lnum2=sum(dnorm(resim4$mu[lulu,],mean=meand,sd=resim4$sig[lulu,],log=TRUE)+
dgamma(resim4$sig[lulu,],10,vard,log=TRUE)-2*log(resim4$sig[lulu,]))+
sum((rep(0.5,mix$k)-1)*log(resim4$p[lulu,]))+
lgamma(sum(rep(0.5,mix$k)))-sum(lgamma(rep(0.5,mix$k)))
lchibapprox4=lnum1+lnum2-log(resim4$deno)
S=readline(prompt="Type <Return> to continue : ")
########## k=5
k=5
mu=rnorm(k,mean=meand,sd=sdd/k)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(1,k))
mix=list(k=k,p=p,mu=mu,sig=sig)
resim5=gibbs(200,datha,mix)
lulu=order(resim5$lolik)[200]
lnum1=resim5$lolik[lulu]
lnum2=sum(dnorm(resim5$mu[lulu,],mean=meand,sd=resim5$sig[lulu,],log=TRUE)+
dgamma(resim5$sig[lulu,],10,vard,log=TRUE)-2*log(resim5$sig[lulu,]))+
sum((rep(0.5,mix$k)-1)*log(resim5$p[lulu,]))+
lgamma(sum(rep(0.5,mix$k)))-sum(lgamma(rep(0.5,mix$k)))
lchibapprox5=lnum1+lnum2-log(resim5$deno)
S=readline(prompt="Type <Return> to continue : ")
########## k=6
k=6
mu=rnorm(k,mean=meand,sd=sdd/k)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(1,k))
mix=list(k=k,p=p,mu=mu,sig=sig)
resim6=gibbs(200,datha,mix)
lulu=order(resim6$lolik)[200]
lnum1=resim6$lolik[lulu]
lnum2=sum(dnorm(resim6$mu[lulu,],mean=meand,sd=resim6$sig[lulu,],log=TRUE)+
dgamma(resim6$sig[lulu,],10,vard,log=TRUE)-2*log(resim6$sig[lulu,]))+
sum((rep(0.5,mix$k)-1)*log(resim6$p[lulu,]))+
lgamma(sum(rep(0.5,mix$k)))-sum(lgamma(rep(0.5,mix$k)))
lchibapprox6=lnum1+lnum2-log(resim6$deno)
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.