[R] Trouble with pmvnorm?

[peter dalgaard]

On 09 Feb 2014, at 10:56 , Paul Parsons <pparsons298 at gmail.com> wrote:

> 
> Many thanks, Peter. Creating a wrapper function for integrand using Vectorize, and then integrating the wrapper, does indeed solve the problem. I tried your final suggestion, but the variable x still gets passed into pmvnorm inside the new mean and variance matrix, leading to the same problem when the integrate function vectorizes x.

You missed the point: There's nothing to integrate if you do it that way. All you need is the marginal distribution of the differences.

-pd

> All the best
> Paul
> 
> On 8 Feb 2014, at 18:04, peter dalgaard wrote:
> 
>> ou 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))??
> 

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