MCMC gradually slows down

[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

______________________________________________
[hidden email] 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.

[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: [phone]
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: [hidden email] [[hidden email]] On Behalf Of Jens Malmros [[hidden email]]
Sent: Sunday, November 08, 2009 8:11 PM
To: [hidden email]
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

______________________________________________
[hidden email] 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.
______________________________________________
[hidden email] 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.

[jholtman]

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 <[hidden email]> 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
>
> ______________________________________________
> [hidden email] 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
[phone]

What is the problem that you are trying to solve?

______________________________________________
[hidden email] 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.

[Jens Malmros]

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 <[hidden email]> 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
>>
>> ______________________________________________
>> [hidden email] 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.
>>
>
>
>

______________________________________________
[hidden email] 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.

[William Dunlap]

> -----Original Message-----
> From: [hidden email]
> [mailto:[hidden email]] On Behalf Of Jens Malmros
> Sent: Sunday, November 08, 2009 11:11 AM
> To: [hidden email]
> 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
>
> ______________________________________________
> [hidden email] 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.
>

______________________________________________
[hidden email] 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.