You almost said it yourself: Your integrand doesn't vectorize. The direct culprit is the following:
If x is a vector, what is lower=c(x,x,x,x)? A vector of length 4*length(x). And pmvnorm doesn't vectorize so it wouldn't help to have lower= as a matrix (e.g., cbind(x,x,x,x)) instead. A straightforward workaround is to Vectorize() your function. Possibly more efficient to put an mapply() of sorts around the pmvnorm call. However, wouldn't it be more obvious to work out the mean and variance matrix of (x1-x5, x2-x5, x3-x5, x4-x5) and then just pmvnorm(... lower=c(0,0,0,0), upper=c(Inf, Inf, Inf, Inf))?? On 06 Feb 2014, at 21:53 , Paul Parsons <pparsons...@gmail.com> wrote: > Hi > > I have a multivariate normal distribution in five variables. The distribution > is specified by a vector of means ('means') and a variance-covariance matrix > ('varcov'), both set up as global variables. > > I'm trying to figure out the probabilities of each random variable being the > smallest. > > So I've made a function: > > integrand<-function(x){ > > #create new mv normal dist, conditional on fixing the > value of element i to x > sig11 <- varcov[-i,-i] > sig12 <- varcov[,i] > sig12 <- sig12[-i] > sig21 <- varcov[i,] > sig21 <- sig21[-i] > sig22 <- varcov[i,i] > mu1 <- means[-i] > mu2 <- means[i] > > muBar <- mu1 + sig12*(x-mu2)/sig22 > sigBar <- sig11 - (sig12) %*% t(sig21)/sig22 > > #now calculate the probability that variable i takes > the value x, > #and that all other variables are bigger than x... > arg <- dnorm(x,means[i],sigma[i]) > arg <- arg*pmvnorm(lower=c(x,x,x,x), > upper=c(10,10,10,10), mean=muBar,sigma=sigBar) > > return(as.numeric(arg)) > } > > Then I need to perform a 1-d integration of this function over all possible > values of x, which gives the total probability of variable i being the > smallest. > > If I use a numerical integration function with explicit looping then this > works fine. But if I try and use a vectorised integrator (such as the > 'integrate' function), to improve performance, then I run into the following > error message: > > Error in checkmvArgs(lower = lower, upper = upper, mean = mean, corr = corr, > : > ‘diag(sigma)’ and ‘lower’ are of different length > > checkmvArgs is a function required by pmvnorm, so I'm fairly sure that's > where the problem lies. diag(sigma) and lower certainly are of the same > length, so not sure at all what's happening here. Has anyone else encountered > this issue? And, if so, do you know the solution? > > Many thanks > Paul > > ______________________________________________ > R-help@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. -- Peter Dalgaard, Professor, Center for Statistics, Copenhagen Business School Solbjerg Plads 3, 2000 Frederiksberg, Denmark Phone: (+45)38153501 Email: pd....@cbs.dk Priv: pda...@gmail.com ______________________________________________ R-help@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.