[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 5
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.
Re: [R] MCMC gradually slows down
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.malm...@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 5
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.
--
Jim Holtman
Cincinnati, OH
+1 513 646 9390
What is the problem that you are trying to solve?
__
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.
Re: [R] MCMC gradually slows down
Try this:
Instead of:
theta = c()
use:
theta - rep(NA, 50)
or however many iterations you want the algorithm to run.
Best,
--
Wolfgang Viechtbauerhttp://www.wvbauer.com/
Department of Methodology and StatisticsTel: +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 5
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.
Re: [R] MCMC gradually slows down
Thanks a lot, this works.
jim holtman wrote:
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.malm...@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 5
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.
