# about random walk Metropolis algorithm

#### zzeta35

##### New Member
Our model is: Y is distributed from a Negative binomial(r,p) where log(p/1-p)=xb.
The prior of b is a multivariate normal(mu0, C0=Diag(k0,k1,...,km)) and prior of r is Geometric(q).

We want to edit a R code of random walk Metropolis algorithm to obtain samples from the joint posterior distribution of (b,r). Use a multivariate Gaussian proposal for b with variance matrix V. and proposal for r is either increase r by 1, decrease r by 1 or leave r unchanged with equal probability.

The problem here is I don't know how to express the proposal for the posterior. I have a similar code which is Poission regression. There is only one parameter beta.

RWM<-function(nits,beta.start,V.prop,lambda.prop,y,X,prior.var) {

p<-length(beta.start)
betas<-matrix(nrow=nits,ncol=p) # to store the MCMC output

n.accept<-0

beta.curr<-beta.start

log.like.curr<-poisson.log.like(beta.curr,X,y) #likelihood function
log.prior.curr<-log.prior(beta.curr,prior.var)

for (i in 1:nits) {
z<-mvn.sim(V.prop)

# print(beta.curr)
beta.prop<-beta.curr+lambda.prop*z # proposed value
# print(beta.prop)

log.like.prop<-poisson.log.like(beta.prop,X,y)
log.prior.prop<-log.prior(beta.prop,prior.var)

# print(log.like.prop)

log.alpha<-log.prior.prop+log.like.prop-log.prior.curr-log.like.curr

u<-runif(1)

if (log(u) <= log.alpha) { # accept
beta.curr<-beta.prop
log.like.curr<-log.like.prop
log.prior.curr<-log.prior.prop
n.accept<-n.accept+1
}

betas[i,]<-beta.curr
}

return(betas)
}

Thank you for any help. This problem afflict me too much.

#### zzeta35

##### New Member
Sorry, it's the first time to post. I put the code in your formal:
Code:
 RWM<-function(nits,beta.start,V.prop,lambda.prop,y,X,prior.var) {

p<-length(beta.start)
betas<-matrix(nrow=nits,ncol=p) # to store the MCMC output

n.accept<-0

beta.curr<-beta.start

log.like.curr<-poisson.log.like(beta.curr,X,y) #likelihood function
log.prior.curr<-log.prior(beta.curr,prior.var)

for (i in 1:nits) {
z<-mvn.sim(V.prop)

# print(beta.curr)
beta.prop<-beta.curr+lambda.prop*z # proposed value
# print(beta.prop)

log.like.prop<-poisson.log.like(beta.prop,X,y)
log.prior.prop<-log.prior(beta.prop,prior.var)

# print(log.like.prop)

log.alpha<-log.prior.prop+log.like.prop-log.prior.curr-log.like.curr

u<-runif(1)

if (log(u) <= log.alpha) { # accept
beta.curr<-beta.prop
log.like.curr<-log.like.prop
log.prior.curr<-log.prior.prop
n.accept<-n.accept+1
}

betas[i,]<-beta.curr
}

return(betas)
}