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.
