Re: [R] MCMC gradually slows down

2009-11-08 Thread 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  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)>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.


[R] MCMC gradually slows down

2009-11-08 Thread 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)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.