[R] MCMC coding problem

[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:

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