R Code 8 / Metropolis Hastings Steps


#sample from a standard normal using MH with a random walk proposals

tuneSD=1.5

M=2000
samp=1:M
counter=0

#starting value
samp[1]=1.2

for (m in 2:M)
{
samp[m]=samp[m-1] #default if not changed below
candidate=rnorm(1,samp[m-1],tuneSD) #note symmetry
ll.cur=dnorm(samp[m-1],0,1,log=T)
ll.cand=dnorm(candidate,0,1,log=T)
prob=min(exp(ll.cand-ll.cur),1)
coin=rbinom(1,1,prob)
if (coin)
{
samp[m]=candidate
counter=counter+1
}
}

accept.rate=sum(counter)/M

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

#How about a three parameter lognormal?
#let's do it together.
set.seed(123)

library('MCMCpack')

I=500
y=.4+rlnorm(I,-.4,sqrt(.2))

#prior
mu.0=0
s2.0=1
a=.1
b=.05

M=50000
psiSD=.01
counter=0

mu=1:M
s2=1:M
psi=1:M

mu[1]=-.3
s2[1]=.3
psi[1]=min(y)-.1

for (m in 2:M)
{
#conditional on psi,
z=log(y-psi[m-1])
#mu
v=1/((I/s2[m-1])+(1/s2.0))
c=(I*mean(z)/s2[m-1])+(mu.0/s2.0)
mu[m]=rnorm(1,v*c,sqrt(v))
#s2
s2.scale=sum((z-mu[m])^2)/2+b
s2[m]=rinvgamma(1,shape=a+I/2,scale=s2.scale)
#your turn, now sample psi given mu and s2
psi[m]=psi[m-1]
cand=rnorm(1,psi[m-1],psiSD)
if(cand>0)
{
full.cur=sum(dlnorm(y-psi[m-1],mu[m],sqrt(s2[m]),log=T))
full.cand=sum(dlnorm(y-cand,mu[m],sqrt(s2[m]),log=T))
prob=min(exp(full.cand-full.cur),1)
if (rbinom(1,1,prob)){psi[m]=cand;counter=counter+1}
}
}

Research Methods 3020