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

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.

First of all, allocate 'theta' to be the final size you need. Every time through your loop you are extending it by one, meaning you are spending a lot of time copying the data each time. Do something like: theta <- numeric(n) and then see how fast it works. On Sun, Nov 8, 2009 at 2:11 PM, Jens Malmros <jens.malmros at gmail.com> wrote:> 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. >-- Jim Holtman Cincinnati, OH +1 513 646 9390 What is the problem that you are trying to solve?

> -----Original Message----- > From: r-help-bounces at r-project.org > [mailto:r-help-bounces at r-project.org] On Behalf Of Jens Malmros > Sent: Sunday, November 08, 2009 11:11 AM > 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()Change that theta <- c() to theta <- numeric(n) and it will go faster. Growing datasets can result in a lot of unnecessary memory management. The original had an obvious quadratic quality in the plot of n vs. MHastings(n,.2,2) and preallocating theta to its final length made it look linear up to n=100000 and at n=100000 the time for the original was 35 seconds vs 3.75 for the preallocated version. Bill Dunlap Spotfire, TIBCO Software wdunlap tibco.com> 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. >