disc-mh

Outline

We will draw from the posterior of a normal mean using a log normal prior.

Sampling model:

\[ r \sim N(\mu,\sigma^2) \]

Prior:

\[ \log(\mu) \sim N(\lambda,w^2) \]

Set values for prior and \(\sigma\):

mbar = .01
lam = log(mbar)
w = 1.3
sigma = .04

Data

rd = read.csv("http://www.rob-mcculloch.org/data/monthly_returns.csv")
n = length(rd[[1]])
rbar = mean(rd[[1]])

Discretization, Prior and Posterior on a grid

## discretize prior
# grid for mu values
mg = seq(from=.001,to=.05,length.out=100)

## prior
priv = dlnorm(mg,lam,w)
priv = priv/sum(priv)

Posterior on grid.

## likelihood and posterior
loglv = (-(mg-rbar)^2) * (n/(2*sigma^2))
logpost = log(priv) + loglv
logpost = logpost -max(logpost)
postDv = exp(logpost)
postDv = postDv/sum(postDv)

Plot prior and posterior:

## plot posterior and prior
rgy = range(c(priv,postDv))
plot(range(mg),rgy,type="n",xlab="mu",ylab="post and pri")
lines(mg,priv,col="red")
lines(mg,postDv,col="blue")

MCMC with independent proposal and explicit transition matrix

Make transition matrix.

ns = length(mg)
P = matrix(0.0,ns,ns) # will be markov transition matrix
Pi = postDv #posterior
for(j in 1:ns) {
   temp = Pi/Pi[j]
   junk = ifelse(temp<1,temp,1)
   P[j,] = junk*(1/ns)
   P[j,j] = 1-sum(P[j,][-j])
}

Draw by iterating \(P\).

#iid
ND = 10000


#MCMC by iterating with P
DrMCMC = rep(0,ND)
ii = 1
DrMCMC[1]=mg[ii]

for(i in 2:ND) {
   ii  = sample(1:ns,size=1,prob=P[ii,])
   DrMCMC[i] = mg[ii]
}

IID draws:

DrIId = sample(mg,size=ND,replace=TRUE,prob=Pi)

Compare draws.

## compare iid with MCMC
qqplot(DrMCMC[100:ND],DrIId[100:ND])
abline(0,1,col="red",lwd=3)


## acf
par(mfrow=c(2,1))
acf(DrMCMC[100:ND])
acf(DrIId[100:ND])

Eigen analysis of P

## check eigen
temp = eigen(t(P))
plot(temp$values)


PPi = temp$vectors[,1]
PPi = PPi/sum(PPi)
plot(Pi,PPi)
abline(0,1,col="red")

MCMC by usual reject/accept method

Compute prior on grid.

getpri = function(mu,mg,priv) {
   ii = which.min(abs(mu-mg))
   return(priv[ii])
}

Iterate chain.

Drv = rep(0,ND)
Drv[1] = mg[1]
for(i in 2:ND) {
   ## proposal
   mup = sample(mg,size=1) #uniform proposal
   ## accept prob
   # log post at proposal
   logpip = (-(mup-rbar)^2) * (n/(2*sigma^2)) + log(getpri(mup,mg,priv))
   # log post at current (i-1)
   logpic = (-(Drv[i-1]-rbar)^2) * (n/(2*sigma^2)) + log(getpri(Drv[i-1],mg,priv))
   pacc = min(1,exp(logpip-logpic))
   ## draw accept
   acc = rbinom(1,1,pacc)
   if(acc) { # if accept
      Drv[i] = mup
   } else { # repeat
      Drv[i] = Drv[i-1]
   }
}