R help - Feb 2014 - Trouble with pmvnorm?

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

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 <pparsons298 at 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 at 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.mes at cbs.dk  Priv: PDalgd at gmail.com