+ Reply to Thread
Results 1 to 4 of 4

Thread: about random walk Metropolis algorithm

  1. #1
    Points: 13, Level: 1
    Level completed: 25%, Points required for next Level: 37

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

    about random walk Metropolis algorithm




    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.

  2. #2
    Devorador de queso
    Points: 95,889, Level: 100
    Level completed: 0%, Points required for next Level: 0
    Awards:
    Posting AwardCommunity AwardDiscussion EnderFrequent Poster
    Dason's Avatar
    Location
    Tampa, FL
    Posts
    12,937
    Thanks
    307
    Thanked 2,630 Times in 2,246 Posts

    Re: about random walk Metropolis algorithm

    I don't have emotions and sometimes that makes me very sad.

  3. #3
    Points: 13, Level: 1
    Level completed: 25%, Points required for next Level: 37

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

    Re: about random walk Metropolis algorithm

    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)
     }

  4. #4
    Points: 13, Level: 1
    Level completed: 25%, Points required for next Level: 37

    Posts
    3
    Thanks
    0
    Thanked 0 Times in 0 Posts

    Re: about random walk Metropolis algorithm


    Does anybody can help?

+ Reply to Thread

           




Tags for this Thread

Posting Permissions

  • You may not post new threads
  • You may not post replies
  • You may not post attachments
  • You may not edit your posts






Advertise on Talk Stats