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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._