[R] MCMC gradually slows down

Viechtbauer Wolfgang (STAT) Wolfgang.Viechtbauer at STAT.unimaas.nl
Sun Nov 8 20:42:35 CET 2009


Try this:

Instead of:

theta = c()

use:

theta <- rep(NA, 500000)

or however many iterations you want the algorithm to run.

Best,

--
Wolfgang Viechtbauer                        http://www.wvbauer.com/
Department of Methodology and Statistics    Tel: +31 (0)43 388-2277
School for Public Health and Primary Care   Office Location:
Maastricht University, P.O. Box 616         Room B2.01 (second floor)
6200 MD Maastricht, The Netherlands         Debyeplein 1 (Randwyck)
________________________________________
From: r-help-bounces at r-project.org [r-help-bounces at r-project.org] On Behalf Of Jens Malmros [jens.malmros at gmail.com]
Sent: Sunday, November 08, 2009 8:11 PM
To: r-help at r-project.org
Subject: [R] MCMC gradually slows down

Hello,

I have written a simple Metropolis-Hastings MCMC algorithm for a
binomial parameter:

MHastings = function(n,p0,d){
        theta = c()
        theta[1] = p0
        t =1
        while(t<=n){
                phi = log(theta[t]/(1-theta[t]))
                phisim = phi + rnorm(1,0,d)
                thetasim = exp(phisim)/(1+exp(phisim))
                r = (thetasim)^4*(1-thetasim)^8/(theta[t]^4*(1-theta[t])^8)
                if(runif(1,0,1)<r){
                        theta[t+1] = thetasim
                } else {
                        theta[t+1] = theta[t]
                }
                t = t+1
                if(t%%1000==0) print(t) # diagnostic
        }
        data.frame(theta)
}

The problem is that it gradually slows down. It is very fast in the
beginning, but slows down and gets very slow as you reach about 50000
iterations and I need do to plenty more.

I know there are more fancy MCMC routines available, but I am really
just interested in this to work.

Thank you for your help,
Jens Malmros

______________________________________________
R-help at r-project.org mailing list
https://stat.ethz.ch/mailman/listinfo/r-help
PLEASE do read the posting guide http://www.R-project.org/posting-guide.html
and provide commented, minimal, self-contained, reproducible code.



More information about the R-help mailing list