# Probabilities: Why does this not yield multinomial logistic regression expression?

Update: The formula

Exp[a]/(Exp[a] + Exp[b] + Exp[c])


results from using ExtremeValueDistribution rather than GumbelDistribution as originally stated.

End of Update

Evaluating the following command:

Probability[(x > y ) && (x > z), {x, y, z} \[Distributed]
ProductDistribution[GumbelDistribution[a, 1],
GumbelDistribution[b, 1], GumbelDistribution[c, 1]]]


yields:

> (E^(2 a) (E^(a + b) + E^(a + c) + 2 E^(b + c)))/((E^a + E^b) (E^a + E^
> c) (E^(a + b) + E^(a + c) + E^(b + c)))


shouldn't it be:

exp[a]/(exp[a] + exp[b] + exp[c])


Binary Logit works fine, however:

Probability[
x > y, {x, y} \[Distributed]
ProductDistribution[GumbelDistribution[a, 1],
GumbelDistribution[b, 1]]]


yielding

E^a/(E^a + E^b)

• Why do you think exp[a]/(exp[a] + exp[b] + exp[c]) is the answer? And by the way, what is MNL and why do you think every potential helper would know what it stands for? Mar 5, 2018 at 17:02
• en.wikipedia.org/wiki/Multinomial_logistic_regression Also, it works fine for the binary logit Probability[ x > y, {x, y} [Distributed] ProductDistribution[GumbelDistribution[a, 1], GumbelDistribution[b, 1]]] Mar 5, 2018 at 17:08

It doesn't yield that formula because it shouldn't. Simulations can provide simple and direct checks:

SeedRandom[12345];
a = 1;
b = 2;
c = 3;
nsim = 1000000;
x = RandomVariate[GumbelDistribution[a, 1], nsim];
y = RandomVariate[GumbelDistribution[b, 1], nsim];
z = RandomVariate[GumbelDistribution[c, 1], nsim];

N[Total[Boole[#[[1]] > #[[2]] && #[[1]] > #[[3]]] & /@
Transpose[{x, y, z}]]/nsim]
(* 0.053482 *)

N[(E^(2 a) (E^(a + b) + E^(a + c) + 2 E^(b + c)))/((E^a + E^b) (E^a + E^c) (E^(a + b) + E^(a + c) + E^(b + c)))]
(* 0.0533853 *)

N[Exp[a]/(Exp[a] + Exp[b] + Exp[c])]
(* 0.0900306 *)