# CHAPTER 6: R COMMANDS
# 22/12/2006

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

# license DATASET IMPLEMENTATION

license=read.table("license")
image(1:25,1:105,as.matrix(license),col=grey(0:255/255))
data=scan("license")
data=jitter(data,10)
data=log((data-min(data)+0.01)/(max(data)+0.01-data))
hist(data,nclass=100,prob=TRUE)

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

# GIBBS SAMPLER FOR A MEAN MIXTURE

# SAMPLE CONSTRUCTION

N=5*100

sampl=c(rnorm(100,sd=.1),rnorm(100,mean=2.9,sd=.1),rnorm(100,mean=-2.9,sd=.1),rnorm(100,mean=-4.5,sd=.1),rnorm(100,mean=4.5,sd=.1))

# GRID

mu1=seq(-5.5,5.5,length=150) 
mu2=mu1 #seq(2.0,3.0,.008)
mo1=mu1%*%t(rep(1,length=length(mu2)))
mo2=(rep(1,length=length(mu2)))%*%t(mu2)
ca1=-0.5*mo1*mo1
ca2=-0.5*mo2*mo2

# LIKELIHOOD SURFACE

like=0*mo1
for (i in 1:N) like=like+log(0.7*exp(ca1+sampl[i]*mo1)+0.3*exp(ca2+sampl[i]*mo2))
like=like+.1*(ca1+ca2)
like=like-min(like)

# GIBBS SAMPLER

Nsim=100000
muz=matrix(0,ncol=2,nrow=Nsim)

gu1=mean(sampl)+rnorm(1)
gu2=mean(sampl)+rnorm(1)
muz[1,]=(c(gu1,gu2))+rnorm(2)

for (t in 2:Nsim){

# ALLOCATION

fact=.3*sqrt(exp(gu1^2-gu2^2))/.7
probs=1/(1+fact*exp(sampl*(gu2-gu1)))
zeds=(runif(N)<probs)

# GIBBS SAMPLING

gu1=rnorm(1)/sqrt(.1+length(zeds[zeds==1]))+(sum(sampl[zeds==1]))/(.1+length(zeds[zeds==1]))
gu2=rnorm(1)/sqrt(.1+length(zeds[zeds==0]))+(sum(sampl[zeds==0]))/(.1+length(zeds[zeds==0]))
muz[t,]=(c(gu1,gu2))

}

# OUTPUT

image(mu1,mu2,like,xlab=expression(mu[1]),ylab=expression(mu[2]),col=heat.colors(250))
contour(mu1,mu2,like,add=T,nlevels=100)
points(muz,pch=18,cex=.5,col="steelblue")

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

# TEMPERING WITH MEAN MIXTURES

# MEAN MIXTURE OF NORMALS

trumyn=2.5
wayt=.8
ssiz=500
sampl=rnorm(ssiz)+(runif(ssiz)<wayt)*trumyn

mu1=seq(min(sampl),max(sampl),.1)
mu2=mu1
mo1=mu1%*%t(mu2/mu2)
mo2=(mu2/mu2)%*%t(mu2)
ca1=-0.5*mo1*mo1
ca2=-0.5*mo2*mo2
like=0*mo1
for (i in 1:ssiz) like=like+log(wayt*exp(ca1+sampl[i]*mo1)+(1-wayt)*exp(ca2+sampl[i]*mo2))
like=like+.1*(ca1+ca2)
like=like-min(like)

# TEMPERING

T=250000

for (alpha in c(1,.3,.1))
{ 	

X11(width=13,height=7)
par(mar=c(4,4,1,1))

tau=.1	

mmu1=rep(1.5,T)
mmu2=rep(3,T) 
olike=sum(log(wayt*exp(-.5*mmu1[1]^2+sampl*mmu1[1])+(1-wayt)*exp(-.5*mmu2[1]^2+sampl*mmu2[1])))

for (t in 2:T)
{

pro1=mmu1[t-1]+tau*rnorm(1)
pro2=mmu2[t-1]+tau*rnorm(1)

newlike=sum(log(wayt*exp(-.5*pro1^2+sampl*pro1)+(1-wayt)*exp(-.5*pro2^2+sampl*pro2)))
if (log(runif(1))/alpha<newlike-olike)
{
mmu1[t]=pro1;mmu2[t]=pro2;olike=newlike
}
else
{
mmu1[t]=mmu1[t-1];mmu2[t]=mmu2[t-1]
}
}

image(mu1,mu2,like,xlab=expression(mu[1]),ylab=expression(mu[2]),col=heat.colors(250),main=paste(alpha))
contour(mu1,mu2,like,levels=seq(min(like),max(like),(max(like)-min(like))/100),add=TRUE)
points(mmu1,mmu2,pch=20,cex=.3)

}

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

# ESTIMATION OF A MIXTURE OF NORMAL DISTRIBUTIONS WITH
# ALL PARAMETERS UNKNOWN

# A MIXTURE WILL BE REPRESNETED AS A LIST WITH COMPONENTS k, mu, sig and p

# THE DATASET IS CALLED datha

datha=data
n=length(datha)

meand=mean(datha)
vard=var(datha)
sdd=sqrt(vard)

rdirichlet=function(n=1,par)
{
k=length(par)
mat=matrix(0,n,k)
for (i in 1:n) 
{
sim=rgamma(k,shape=par,scale=1)
mat[i,]=sim/sum(sim)
}
mat
}

# INITIALISATION

k=3
mu=rnorm(k,mean=meand,sd=sdd/3)
sig=1/rgamma(k,shape=10,scale=vard)
p=rdirichlet(par=rep(5.5,k))
mix=list(k=k,p=p,mu=mu,sig=sig)

# DIRICHLET LOG-PRIOR (1/2,....,1/2)

Dirichlet=function(p)
{
k=length(p)
ou=-0.5*sum(log(p))+lgamma(k*.5)-k*lgamma(.5)
ou
}

# OVERALL LOG-PRIOR

prior=function(mix,meanp=meand,sdp=sdd,varp=vard)
{
ou=Dirichlet(mix$p)
for (i in 1:mix$k) ou=ou+dnorm(mix$mu[i],mean=meanp,sd=sdp,log=T)+dgamma(1/mix$sig[i],shape=3,rate=varp,log=T)
ou
}

# LOG-LIKELIHOOD

like=function(mix)
{
ou=0
if (n>0)
{
for (i in 1:n) ou=ou+log(sum(mix$p*dnorm(datha[i],mean=mix$mu,sd=sqrt(mix$sig))))
}
ou
}

# LOG-POSTERIOR

lpost1=function(mix,meanp=meand,sigp=sdd,varp=vard)
{
like(mix)+prior(mix,meanp,sigp,varp)
}

# GIBBS SAMPLING

gibbs=function(niter,max,meanp=meand,sigp=sdd,varp=vard)
{

mix=max
z=rep(0,mix$k)
ssiz=rep(0,mix$k)
nxx=ssiz
ssum=ssiz
mug=matrix(0,nrow=niter,ncol=mix$k)
sigg=mug
prog=mug

for (i in 1:niter)
{
for (t in 1:n)
{
prob=mix$p*dnorm(datha[t],mean=mix$mu,sd=sqrt(mix$sig))
z[t]=sample(x=1:mix$k,size=1,prob=prob)
}
for (j in 1:mix$k)
{
ssiz[j]=sum(z==j)
nxj[j]=sum(as.numeric(z==j)*datha)
}
mug[i,]=rnorm(mix$k,mean=((3*meanp/varp)+nxj/mix$sig)/((3/varp)+(ssiz/mix$sig)),sd=1/sqrt((3/varp)+(ssiz/mix$sig)))
mix$mu=mug[i,]
for (j in 1:mix$k) ssum[j]=sum(as.numeric(z==j)*(datha-mix$mu[j])^2)
sigg[i,]=1/rgamma(mix$k,shape=10+0.5*ssiz,rate=varp+0.5*(mix$mu-meanp)^2+0.5*ssum)
mix$sig=sigg[i,]

prog[i,]=rdirichlet(1,par=ssiz+0.5)
mix$p=prog[i,]

}
list(k=k,mu=mug,sig=sigg,p=prog)
}

# BAYESIAN ESTIMATION OF A MIXTURE OF NORMAL DISTRIBUTIONS
# WITH AN UNKNOWN NUMBER OF COMPONENTS
# POISSON P(4) PRIOR ON k

Kmax=15
lambda=4

# SINGLE SIMULATION LOG-PRIOR

singleprior=function(mix,i)
{
ou=dbeta(mix$p[i],.5,.5*(mix$k-1),log=T)+
dnorm(mix$mu[i],mean=meand,sd=sdd,log=T)+
dgamma(1/mix$sig[i],shape=3,rate=vard,log=T)
ou
}

# BASIC BIRTH-AND-DEATH MOVE

BD=function(mix)
{
if ((mix$k<Kmax)&&(mix$k>1)) { kprop=mix$k+sample(c(-1,1),1) }
else
{
if (mix$k==Kmax) { kprop=(mix$k-1) }
else
{
kprop=2
}
}
if (kprop==(mix$k-1))
{ 
dead=sample(1:mix$k,1) 
propmix=list(k=kprop,mu=mix$mu[-dead],sig=mix$sig[-dead],p=mix$p[-dead]/(1-mix$p[dead]))
jacob=-(kprop-1)*log(1-mix$p[dead]) + singleprior(mix,dead)
}
if (kprop==(mix$k+1))
{
propp=rbeta(1,1,kprop)
propp=c((1-propp)*mix$p,propp)
propmix=list(k=kprop,mu=c(mix$mu,rnorm(1,mean=meand,sd=sdd)),sig=c(mix$sig,1/rgamma(1,shape=3,rate=vard)),p=propp)
jacob=(kprop-2)*log(1-propp[kprop]) - singleprior(propmix,kprop)
  }
propacpt=dpois(kprop,lambda,log=T)-dpois(mix$k,lambda,log=T)+jacob+
log(1+(kprop==1)+(kprop==Kmax))+prior(propmix)+like(propmix)-
log(1+(mix$k==1)+(mix$k==Kmax))-prior(mix)-like(mix)

if (log(runif(1))<propacpt) mix=propmix

mix
}

# THE RJMCMC

RJMCMC=function(mix,T)
{
kaze=rep(0,1000)
mus=mix$mu
sigs=mix$sig
pees=mix$p

for (t in 1:1000)
{
for (ins in 1:(T/1000)) mix=BD(mix)
kaze[t]=mix$k
mus=c(mus,mix$mu)
sigs=c(sigs,mix$sig)
pees=c(pees,mix$p)
}
list(k=kaze,mu=mus,sig=sigs,p=pees,last=mix)
}

# DRAWING THE MIXTURE

mox=function(mix)
{
hist(datha,prob=T,main="",xlab="",nclass=100,col="gold")
x=seq(min(datha),max(datha),length=1000)
y=mix$p[1]*dnorm(x,mean=mix$mu[1],sd=sqrt(mix$sig[1]))
for (i in 2:mix$k) y=y+mix$p[i]*dnorm(x,mean=mix$mu[i],sd=sqrt(mix$sig[i]))
lines(x,y,lwd=2,col="sienna3")
for (i in 2:mix$k) lines(x,mix$p[i]*dnorm(x,mean=mix$mu[i],sd=sqrt(mix$sig[i])),col="chocolate",lty=3)
}

# ALLOCATING THE OBSERVATIONS

alloc=function(mix)
{
missing=rep(0,n)
for (i in 1:n)
{
prob=mix$p*dnorm(datha[i],mean=mix$mu,sd=sqrt(mix$sig))
missing[i]=order(prob)[mix$k] 
}
missing
}