> -----Original Message-----
> From: r-help-boun...@r-project.org 
> [mailto:r-help-boun...@r-project.org] On Behalf Of Jens Malmros
> Sent: Sunday, November 08, 2009 11:11 AM
> 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()

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

Reply via email to