[R] pmnorm: probabilites don't sum up to 1
nicolas.berkowitsch at unibas.ch
Wed Dec 15 16:50:48 CET 2010
Thx for your comment - I guess you're right... I made a drawing to
better understand which quadrants overlap
So I ende up with this Code:
library(mnormt) # can handle multivariate normal distributions
uX = 2
uY = 1
uZ = .5
mean = c(uX, uY, uZ)
LX = matrix(c(1,-1,0,1,0,-1), 2, 3, byrow = TRUE)
LY = matrix(c(-1,1,0,0,1,-1), 2, 3, byrow = TRUE)
LZ = matrix(c(-1,0,1,0,-1,1), 2, 3, byrow = TRUE)
muX = LX %*% mean
muY = LY %*% mean
muZ = LZ %*% mean
Sigma = diag(2)
mean = c(0,0)
pX = pmnorm(LXmean,c(0,0),varcov)/(1-(pmnorm(-LXmean,c(0,0),varcov)))
pY = pmnorm(LYmean,c(0,0),varcov)/(1-(pmnorm(-LYmean,c(0,0),varcov)))
pZ = pmnorm(LZmean,c(0,0),varcov)/(1-(pmnorm(-LZmean,c(0,0),varcov)))
pX + pY + pZ
pX is now expressed as a relation:
uX-uY>0 & uX-uZ> = pmnorm(LXmean,c(0,0),varcov) --> counter
in relation to correspondig surface area =
(1-(pmnorm(-LXmean,c(0,0),varcov))) --> denominator
This now sums up to 1 (or nearly) - I guess this should be the correct
way of doing it, right?
Zitat von peter dalgaard <pdalgd at gmail.com>:
> On Dec 15, 2010, at 11:40 , Nicolas Berkowitsch wrote:
>> Dear list member,
>> I struggle with the problem, why the probabilities of choosing one
>> of three mutually exclusive alternatives don?t sum up to 1!
>> Let?s assume we have three alternatives X, Y, and Z. Let?s further
>> assume we know their respective utilities:
>> uX, uY, uZ. I?m interested in calculating the probability of
>> choosing X, Y, and Z.
>> Since I assume that the alternatives are mutually exclusive, the
>> probabilities p(X), p(Y), and p(Z) have to sum up to one. The
>> utilities of the 3 alternatives can be expressed in 2 utility
>> differences and, hence, the multivariate case reduces to a
>> bivariate normal distribution. If I assume that X, Y, and Z are
>> independent, their corresponding correlations have to be zero and,
>> hence, the variance-covariance-matrices are set to be a
>> diagonal-matrix (i.e., identity-matrix).
>> To calculate p(X), p(Y), and p(Z) I was using the following R-code:
>> library(mnormt) # can handle multivariate normal distributions
>> uX = 2
>> uY = 1
>> uZ = .5
>> mu = c(uX, uY, uZ)
>> LX = matrix(c(1,-1,0,1,0,-1), 2, 3, byrow = TRUE)
>> LY = matrix(c(-1,1,0,0,1,-1), 2, 3, byrow = TRUE)
>> LZ = matrix(c(-1,0,1,0,-1,1), 2, 3, byrow = TRUE)
>> muX = LX %*% mu
>> muY = LY %*% mu
>> muZ = LZ %*% mu
>> Sigma = diag(2)
>> mean = c(0,0)
>> pX = pmnorm(muX, mean, Sigma)
>> pY = pmnorm(muY, mean, Sigma)
>> pZ = pmnorm(muZ, mean, Sigma)
>> pX + pY + pZ
>> I don?t see why the three probabilities don?t sum up to 1?
>> I know two ?solutions? to this problem so far. However, neither of
>> them satisfies me:
>> 1. I can set pZ to 1 ? pX ? pY, but doing so, returns a different
>> result for pZ, than calculating pZ directly using pmnorm.
>> 2. I could calculate the relationship of pX to the sum of pX + pY +
>> pZ (? pX/(pX + pY + pZ) )
>> Can anyone explain to me why the probabilities don?t sum up to 1?
>> How should I rewrite the R-code to overcome this problem?
>> Thanks a lot for any advice!
> I don't think the pX, pY, pZ are probabilities of choosing X, Y, Z.
> If you think that they are, then you need to explain it more
> What they are are probabilities of three lower-left quadrants with
> different origins (muX, muY, muZ). Such quadrants will in general
> overlap, so there is no reason to expect their probabilities to sum
> to any particular value. If you were expecting to have a partition
> of 2d space into three disjoint regions and calculating the
> probability of each region, then pmnorm is not the right tool.
> Peter Dalgaard
> 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
lic. phil. Nicolas A. J. Berkowitsch
Fakultät für Psychologie
Tel. +41 61 267 05 75
E-Mail nicolas.berkowitsch at unibas.ch
This message was sent using IMP, the Internet Messaging Program.
More information about the R-help