Skip to content

MCMC gradually slows down

5 messages · Jens Malmros, Viechtbauer Wolfgang (STAT), jim holtman +1 more

Original
Jens Malmros
# 
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
Viechtbauer Wolfgang (STAT)
# 
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)
jim holtman
# 
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:

  
    

      
      
    

  


  


      

    

    

      
      

        


  

        

      Jens Malmros
    

    

      
      #
    

  


  

      
Thanks a lot, this works.

        
jim holtman wrote:

            

      
      
      
    

        

      
      
    

            

      
      
      
    

  
  


  


      

    

    

      
      

        


  

        

      William Dunlap
    

    

      
      #
    

  


  

          

      
      
      
    

        
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