[R] MCMC coding problem

[ben@zoo.ufl.edu]

   I used the traceback() command (run it immediately after you encounter
an error in R) to see what was happening.

The result was

3: rep(data, t1)
2: array(0, dim = c(3, d, len))
1: mcmcmd(2, 0.6, 2, 2, 0.3, 100)

which says that the problem is in the "array()" command in the 6th line of
your code.  Looking more closely, it seems that the problem is that len
is not a positive integer unless ntimes is a multiple of 1000 greater than
1000.  Perhaps S-PLUS automatically rounds the value to something
sensible?  (Giving an error instead is in keeping with the R philosophy of
giving an error when something weird happens, rather than covering it up
and keeping going.)

(When I try mcmcmd(2,0.6,2,2,0.3,2000), I get a bit farther, to

Error: couldn't find function "rmvnorm"

which is probably a function you have defined elsewhere.)

 good luck,
   Ben Bolker

On Thu, 30 Aug 2001, Monnie McGee wrote:

> Dear All,
>
> I am trying to convert some S-plus code that I have to run MCMC into
> R-code.  The program works in S-plus, but runs slowly.
>
> I have managed to source the program into R.  R recognizes that the program
> is there; for example, it will display the code when I type the function
> name at the prompt.  However, the program will not run.  When I try to run
> the program, I get the following error message:
>
> > error in rep(data,t1): invalid number of copies in "rep".
>
> This happens no matter what values that I give the parameters when I call
> the function.
>
> Any suggestions on clearing this up would be appreciated!  The function is
> printed below.
>
> Many Thanks,
> Monnie McGee
>
> mcmcmd <- function(d, w, mean, phi2, varwt, ntimes)
> {
> # A function to run a multi dimensional MCMC for a bivariate
> # normal distribution
> # d = dimension of problem
> # w = weight assigned to larger of two modes in target distribution
> # mu = mean of second mode
> # phi2 = multiplier of variance factor 2.38
> # varwt = weight (0 < varweight <= 1) applied to variance of second part of
> mixture.
> # ntimes = number of iterations
> #
> # phi is the standard deviation of the proposal distribution.
> # y is the data matrix, ystar is the proposal state
>
> phi <- (2.38^2)/d * diag(phi2,d)
> y <- matrix(0,ncol=d,nrow=ntimes)
> ystar <- matrix(0,ncol=d,nrow=1)
> len <- ntimes/1000 - 1
> meanmat <- matrix(0,nrow=len,ncol=d)
> quantmat <- array(0, dim = c(3,d,len))
> count <- 0
>
> # print("Comment out the browser before running with ntimes large!!")
>
> browser()
> for (i in 1:(ntimes-1))
> 	{
> 	ystar <- rmvnorm(n=1,mu=y[i,],sigma=phi)
> 	numer <- (w*dmvnorm(ystar,mu=rep(0,d),sigma=diag(0,d)) +
> (1-w)*dmvnorm(ystar,mu=rep(mean,d),sigma=diag(varwt,d)))*dmvnorm(y[i,],mu=ys
> tar,sigma=phi)
> 	denom <- (w*dmvnorm(y[i,],mu=rep(0,d),sigma=diag(0,d)) +
> (1-w)*dmvnorm(y[i,],mu=rep(mean,d),sigma=diag(varwt,d)))*dmvnorm(ystar,mu=y[
> i,],sigma=phi)
> 	prob <- numer/denom
> 	alpha <- min(prob,1)
> 	u <- runif(1)
> 	#browser()
> 	if (u <= alpha) {y[i+1,] <- ystar;count <- count+1} else y[i+1,]<- y[i,]
> 	j <- i/1000
> 	if (j == ceiling(j))
> 	{
> 	quantmat[,,j] <- apply(y,2,quantile,probs=c(.025,.5,.975))
> 	meanmat[j,] <- apply(y,2,mean)
> 	}
> 	}
> return(y,count,quantmat,meanmat)
> }
>
>
> -.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
> r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
> Send "info", "help", or "[un]subscribe"
> (in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
> _._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
>

-- 
318 Carr Hall                                bolker at zoo.ufl.edu
Zoology Department, University of Florida    http://www.zoo.ufl.edu/bolker
Box 118525                                   (ph)  352-392-5697
Gainesville, FL 32611-8525                   (fax) 352-392-3704

-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-.-
r-help mailing list -- Read http://www.ci.tuwien.ac.at/~hornik/R/R-FAQ.html
Send "info", "help", or "[un]subscribe"
(in the "body", not the subject !)  To: r-help-request at stat.math.ethz.ch
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._