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-boun...@r-project.org [r-help-boun...@r-project.org] On Behalf Of Jens Malmros [jens.malm...@gmail.com] Sent: Sunday, November 08, 2009 8:11 PM To: r-help@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@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. ______________________________________________ R-help@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.