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. >