[R] Gibbs sampler...did it work?

jim holtman jholtman at gmail.com
Mon Jan 26 00:57:01 CET 2009


It is not that you are out of memory; one of your two 's2yg' objects
is not large enough (improperly dimensioned?) so you get the subscript
error.  You can put the following in your script to catch the error
and then examine the values:

options(error=utils::recover)


On Sun, Jan 25, 2009 at 6:25 PM, ekwaters <ekwaters at unimelb.edu.au> wrote:
>
> I am writing a Gibbs sampler. I think it is outputting some of what I want,
> in that I am getting vector of several thousand values (but not 10,000) in a
> txt file at the end.
>
> My question is, is the error message (see below) telling me that it can't
> output 10,000 values (draws) because of a limitation in my memory, file
> size, shape etc, or that there is an error in the sampler itself?
>
>> s2eg2=1/rgamma(mg2,(12/2),.5*t(residuals(lm(yg[,1]~xg-1))%*%residuals(lm(yg[,1]~xg-1))))
>> for(i in 1:mg2){
> + s2yg[i,]=parsy+t(rnorm(1,mean=0,sd=s2ygscale[i])%*%chol(s2eg2[i]*xgtxgi))
> + write(c(s2yg[i,],s2eg2[i]),
> + file="/media/DataTravelerMini/KINGSTON/Honours/R/IPR/s2yg2.txt", append=T,
> ncolumns=1)
> + if(i%%50==0){print(c(s2yg[i,],s2eg2[i]))}}
>
> I GET A BUNCH OF NUMBERS PRINTED HERE, THE OUTPUTTED VALUES WHICH ALSO
> APPEAR IN A TEXT FILE.  I HIT ABOUT 2000 VALUES, THEN I GET THIS MESSAGE:
>
> Error in s2yg[i, ] = parsy + t(rnorm(1, mean = 0, sd = s2ygscale[i]) %*%  :
>        subscript out of bounds
>
>
> Ned
>
>
> --
> View this message in context: http://www.nabble.com/Gibbs-sampler...did-it-work--tp21658246p21658246.html
> Sent from the R help mailing list archive at Nabble.com.
>
> ______________________________________________
> R-help at 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?




More information about the R-help mailing list