# Calculate the mean of ordered distribution of the defined distribution

I am trying to evaluate the mean of the following:

Let $a\sim Binom(10,1/2)$ and consider the random variable

$\quad \quad d=\left|\frac{2}{\sqrt{10}}a-\sqrt{10}\right|$

Let $d_1,\dots,d_{20}$ be i.i.d copies of $d$, and $d_{(1)},\dots, d_{(20)}$ be the non-decreasing rearrangements, I want to calculate the mean of $d_{(k)}$, say for example, $k=15$.

I used the following codes:

d = TransformedDistribution[Abs[Sqrt[2/5]a-Sqrt[10]],
a\[Distributed]BinomialDistribution[10,1/2]]

f = OrderDistribution[{d, 20}, 15]

Mean[f]


But the output was

Mean[OrderDistribution[{TransformedDistribution[
Abs[-Sqrt[10] + Sqrt[2/5] \[FormalX]], \[FormalX] \[Distributed]
BinomialDistribution[10, 1/2]], 20}, 15]]


No error, no answer. I am confused. Anybody know the correct way to do this? I 'll be really appreciated for the help.

• There's a syntax error in your definition of d (spurious brackets around the second Sqr). – Sjoerd C. de Vries Sep 3 '15 at 21:42
• @SjoerdC.deVries I see. I fixed that, but still didn't work – Frank Lu Sep 3 '15 at 21:47
• I suspect that the TransformedDistribution expression is too complex for Mathematica. You might try using just a as an argument to see whether at least your syntax is correct. My guess would be the Abs is a problem. – Sjoerd C. de Vries Sep 3 '15 at 21:52
• Ok, there's a workaround. Presented here for a slightly simpler case: d = TransformedDistribution[Abs[a - 1], a \[Distributed] BinomialDistribution[10, 1/2]]; newD = ProbabilityDistribution[PDF[d, x], {x, 0, 10, 1}]; Mean@OrderDistribution[{newD, 20}, 15] – Stefan R Sep 3 '15 at 22:16

I'm a tad surprised MMA has difficulty here, but easy enough to derive:

Clear[z, \[ScriptF], \[ScriptCapitalF], n, j, do]

(* the distribution at hand*)
do = TransformedDistribution[Abs[Sqrt[2/5] a - Sqrt[10]],
a \[Distributed] BinomialDistribution[10, 1/2]];

(* get PDF, pull parts, create pdf/cdf lists*)
pdf = PDF[do, z];
pdf = SortBy[Transpose@{pdf[[1, All, 2, 2]], pdf[[1, All, 1]]}, First];
cdf = Transpose[{pdf[[All, 1]], Accumulate@pdf[[All, 2]]}];

(* make pmf/cmf functions *)
(\[ScriptF][#[[1]]] = #[[2]]) & /@ pdf;
(\[ScriptCapitalF][#[[1]]] = #[[2]]) & /@ cdf;
\[ScriptCapitalF][_] = \[ScriptF][_] = 0;

(* order statistic fn - given n samps from X, kth order, and x *)
(* returns PMF for X=x *)
p[n_, k_, x_] :=
Block[{Power}, Power[_, 0] = 1;
Sum[Binomial[n, j] ((1 - \[ScriptCapitalF][x])^j (\[ScriptCapitalF][x])^(n - j)
- (1 - \[ScriptCapitalF][x] + \[ScriptF][x])^j
(\[ScriptCapitalF][x] - \[ScriptF][x])^(n - j)),
{j, 0,  n - k}]];


And a quick check:

(* test - pull probs and valid values for sim, do table of all orders *)
(* using function and sim... say for 20 RV from dist *)

{ps, vals} = {pdf[[All, 2]], pdf[[All, 1]]};

TableForm[
Table[{n, (p[20, n, #] & /@ vals)*vals // Tr //
N, (Sort /@ RandomChoice[ps -> vals, {100000, 20}])[[All, n]] //
Mean // N}, {n, 20}], TableSpacing -> {1, 2},
TableHeadings -> {None, {"Order", "Fn", "Sim"}}]


Validate our order statistic fn is a proper PMF:

PDF[do, vals] == (p[1, 1, #] & /@ vals)

(* True *)

DeleteDuplicates@Table[(p[20, n, #] & /@ vals) // Tr, {n, 20}]

(* {1} *)


Looks good... sew it up as a function that takes distribution, count, etc., good to go.

BTW, for your k of 15 case, a one-shot use:

(p[20, 15, #] & /@ vals)*vals // Tr


(67159300657458079766884744755289027591228019030511847 Sqrt[5/2]) /95780971304118053647396689196894323976171195136475136

N.b.: You'll likely need to muck with the pmf/cmf parts-puller for different distributions to pull values and corresponding probabilities when there are probabilities with multiple values... if/when time permits, I'll wrap this all up into a single function that does all the work.

Edit: I'd consider reporting this behavior to WRI, IMO it should be able to handle it, but looking at some intermediate results under the covers, something funny seems to be going on.

That said, with a little background there's no need to resort to pulling out the wallet, you can in fact get results for arbitrary n and k on your distribution with not too much effort. Based on the same information about the domain derived above...

d = TransformedDistribution[Abs[Sqrt[2/5] a - Sqrt[10]],
a \[Distributed] BinomialDistribution[10, 1/2]];
pdo = OrderDistribution[{d, n}, k];
pdfo[m_] = CDF[pdo, m] - CDF[pdo, m - Sqrt[2/5]];

mean = Tr@Table[pdfo[m]*m, {m, 0, Sqrt[10], Sqrt[2/5]}];

mean /. {n -> 20, k -> 15}


(67159300657458079766884744755289027591228019030511847 Sqrt[5/2]) /95780971304118053647396689196894323976171195136475136

Matching, of course, the earlier result, and allowing arbitrary n and k. You can of course take mean and do a FullSimpify etc. to get the analytic form...

• Thanks a lot! This solution works perfectly well. Awesome!. – Frank Lu Sep 4 '15 at 19:06
• @FrankLu glad to help. I'll update when I can with a self-contained all-in-one function. And do get your free copy of mathstatica - it has some interesting functionality and a nice little book - I no longer use it (does nothing mma can't do with a little code) but saves time for some problems when one does not have background/time to code same. – ciao Sep 4 '15 at 19:50
• Thanks for the kind suggestion. As I am a beginner of programming, it might take some time for me to get used of theses stuff. I am also curious whether there's a program which works for arbitrary distribution and arbitrary large amount of data, and that would greatly save time. – Frank Lu Sep 4 '15 at 20:19
• Re Edit ---> much improved!! This is, I think, in itself a beautiful illustration of the advantage of openness and friendly interaction with others. All the better for the ability to compare different output and different ways of working, which is what this site is all about. – wolfies Sep 6 '15 at 15:22

OP asks: calculate the mean of $d_{(k)}$ (i.e. the mean of the $k^{th}$ order statistic)

ciao (aka rasher) has provided a nice attempt at a manual workaround. Unfortunately, the solution he provides does not actually work for unknown $k$ or $n$. It gets stuck at the order statistic stage. For example, his function p[n,k,x] returns:

The OP asks about a general solution for the mean of the $k^{th}$ order statistic. One can, in fact, obtain an exact symbolic solution (as a function of $k$ and $n$) ...

The Problem

Given: $\quad \quad \quad A \sim Binomial(10,\frac12) \quad$ and $\quad X=\left|\frac{2}{\sqrt{10}}A-\sqrt{10}\right|$

As the OP notes, the pmf of $X$ can be found with:

PDF[TransformedDistribution[Abs[Sqrt[2/5]a-Sqrt[10]],
a \[Distributed] BinomialDistribution[10, 1/2]], x]


The above is just a list of 6 $x$-values that define the domain of support:

     xvals = Range[0, 5] Sqrt[2/5]


... and a corresponding list of 6 probabilities $f(x) = P(X=x)$, which we can summarise in a List Form, as:

         f  = {126, 210, 120, 45, 10, 1}/512;
domain[f] = {x, xvals}   &&  {Discrete};


We now have a List Form expression of the pmf of $X$.

The Order Statistic Problem

Let $(X_1, \dots, X_n)$ denote a sample of size $n$ drawn on $X$, with parent pmf $f(x)$. Then, the pmf of the $k^{th}$ order statistic, say $g(x_{(k)}) = P(X_{(k)}=x_{(k)})$, corresponding to the same 6 $x$-values that define the domain of support, is:

where I am using the OrderStat function from the mathStatica package for Mathematica. The domain of support is necessarily the same as the parent random variable $X$, i.e. our same 6 $x$-values:

The mean of the $k^{th}$ order statistic, namely $E[X_{(k)}]$ is then:

All done.

• And you'll be sending the OP his required software used here few of charge yes? I for one would appreciate a notation that this is, unless I am mistaken, software you are involved in, and profit from. – ciao Sep 4 '15 at 14:46
• Sure - if he would like a copy, I'd be happy to send him one. The notation you seek is very clearly stated on my user page: don't think you can miss it there :) – wolfies Sep 4 '15 at 14:53