Of course, sorry... The approach below is indeed a much slicker way to
calculate the probability of the ith variable being the smallest:
#construct a new mv norm distribution, for y = x - x_i
#define y = Bx, where B is a matrix
B <- matrix(rep(0,25),nrow=5)
diag(B) <- rep(1,5)
B[,i] <- -1
B <- B[-i,]
newMeans <- as.vector(B %*% means)
newVarcov <- B %*% varcov %*% t(B)
#probability of component i being the smallest...
prob[i] <- as.numeric(pmvnorm(lower=c(0,0,0,0), upper=c(Inf, Inf, Inf,
Inf), mean=newMeans ,sigma=newVarcov ))
Thanks again!
Paul
On Mon, Feb 10, 2014 at 6:42 PM, peter dalgaard <pda...@gmail.com> wrote:
>
> On 09 Feb 2014, at 10:56 , Paul Parsons <pparsons...@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....@cbs.dk Priv: pda...@gmail.com
>
>
>
>
>
>
>
>
>
[[alternative HTML version deleted]]
______________________________________________
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.
