# How can I calculate the probability of a condition concerning ranked values from the same distribution?

How can I calculate in Mathematica the probability of a condition concerning ranked values from the same distribution?

For instance: Pr[x1 - x2 < 1 | {x1, x2} ~ N(0,1)], where x1 is the second ranked (highest value) and x2 the first ranked (lowest value) both drawn from a standard normal distribution.

I'm actually looking at a much more complicated condition* than this, but a general solution applied to the example above should already help me!

*Pr[x1 - x2 > (2 a y + 3 (x2 - L))/(4 a b - 3) | {x1, x2} ~ N(L,var_x), y ~ N(0,var_y)], where x1 is the second ranked and x2 the first ranked drawn from the same distribution.

• Have you already seen OrderDistribution[]? Jan 6, 2017 at 14:19
• I didn't! Thanks. I now tried this, but it keeps giving a probability equal to 1. Any idea what's going wrong? Probability[ x1 - x2 < 1, {x1, x2} [Distributed] OrderDistribution[{NormalDistribution[0, 1], 2}, {1, 2}]] -- (Flipping the k^th order statistic such that x2 is the first ranked (lowest value) and x1 is the second ranked (highest value) does not seem to solve the problem.) Jan 6, 2017 at 14:29
• Shouldn't that be x2 - x1 < 1, considering the ranking of the result of OrderDistribution[]? It is possible that Mathematica does not know a closed form for the underlying integral; you will then need to use NProbability[] instead. Jan 6, 2017 at 14:33
• You're right on the ranking. It should read 'Probability[ x1 - x2 < 1, {x1, x2} [Distributed] OrderDistribution[{NormalDistribution[0, 1], 2}, {2, 1}]]'. Unfortunately this also does not work (gives 0). Problem remains with NProbability (once specifying an accuracy goal). Jan 6, 2017 at 14:59
• What does NProbability[x2 - x1 < 1, {x1, x2} \[Distributed] OrderDistribution[{NormalDistribution[], 2}, {1, 2}]] return for you? Jan 6, 2017 at 15:00

One can take this example further by working outside of the limitations of black boxes, by working with the actual input and outputs involved.

Short version: The exact answer is: Erf[1/2]

Longer version: Given parent random variable $$X \sim N(0,1)$$ with pdf $$f(x)$$:

(source: tri.org.au)

Then, given a sample of size 2 drawn on parent $$X$$, the joint pdf of the $$1^\text{st}$$ and $$2^\text{nd}$$ order statistics is say $$g(x_1,x_2)$$:

(source: tri.org.au)

where I am using the OrderStat function from the mathStatica package for Mathematica here which is the most familiar to me as one of the authors, or one can use Mathematica's in-built OrderDistribution function together with PDF. Either way, viewing output and functional forms is crucial.

Mathematica's in-built Probability function returns Mathematica's Probability function (the same user input), so I will use mathStatica's Prob function here. In particular, we seek $$P(X_2 - X_1 <1)$$:

(source: tri.org.au)

Mathematica does not return a closed form solution, but by being able to see the output, we can take matters further. The Erf[x] term integrates out to 0, leaving an expression equivalent to:

$$E\big[ \text{Erf}\big[\frac{1+X}{\sqrt{2}}\big] \big] \; = \; E\big[ 2 F(1+X) \big] - 1$$

where $$F$$ denotes the cdf of a standard Normal (and where the expectation operator is with reference to parent pdf $$f(x)$$). The problem has now reduced to finding the expectation of the cdf of a standard Normal. Mathematica seems to have trouble with this too, but the expectation does have a closed form ... see, for instance, https://math.stackexchange.com/questions/449875/expected-value-of-normal-cdf ... from which it follows that the closed-form solution is:

$$2 F\big( \frac{1}{\sqrt{2}} \big) - 1 \quad = \quad \text{Erf}[\frac12]$$

which is approximately 0.5205.