Skip to content

MCMC coding problem

2 messages · Monnie McGee, Ben Bolker

Original
Monnie McGee
# 
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:
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
_._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._._
Ben Bolker
# 
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: