5
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ There's a syntax error in your definition of d (spurious brackets around the second Sqr). $\endgroup$ – Sjoerd C. de Vries Sep 3 '15 at 21:42
  • $\begingroup$ @SjoerdC.deVries I see. I fixed that, but still didn't work $\endgroup$ – Frank Lu Sep 3 '15 at 21:47
  • $\begingroup$ 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. $\endgroup$ – Sjoerd C. de Vries Sep 3 '15 at 21:52
  • 2
    $\begingroup$ 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] $\endgroup$ – Stefan R Sep 3 '15 at 22:16
7
$\begingroup$

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"}}]

enter image description here

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...

$\endgroup$
  • $\begingroup$ Thanks a lot! This solution works perfectly well. Awesome!. $\endgroup$ – Frank Lu Sep 4 '15 at 19:06
  • 1
    $\begingroup$ @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. $\endgroup$ – ciao Sep 4 '15 at 19:50
  • 1
    $\begingroup$ 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. $\endgroup$ – Frank Lu Sep 4 '15 at 20:19
  • $\begingroup$ 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. $\endgroup$ – wolfies Sep 6 '15 at 15:22
3
$\begingroup$

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:

enter image description here

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]

enter image description here

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:

enter image description here

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:

enter image description here

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

enter image description here

All done.

$\endgroup$
  • 1
    $\begingroup$ 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. $\endgroup$ – ciao Sep 4 '15 at 14:46
  • $\begingroup$ 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 :) $\endgroup$ – wolfies Sep 4 '15 at 14:53
  • 1
    $\begingroup$ While that's kind of you, I'm sure most users don't check profiles of answerers to determine financial interests, nor should they have to. IMHO, advertisements by example are still... advertisements. This is a site about Mathematica, it's not called nonfreeMathematicAddons.stackechange for a reason. Barring the OPs of any question answered with non-free add-on software getting it free of charge and/or clearly stating it is such in the answer, I find the practice distasteful and not in the spirit of the site. I believe I've made my point. $\endgroup$ – ciao Sep 4 '15 at 15:16
  • 2
    $\begingroup$ But this is a slippery slide you have embarked upon. Should we require all those who work for Wolfram (and who we are very lucky to have here as contributors) to post a disclaimer on their SE posts, because they have an interest in, and/or may derive financial benefit from Mathematica?? It is also a tricky issue to nitpick about, since the biggest gripe about Mathematica (esp from non-Mma users) is the fact that it is not free. Anyway- this is all off-topic, and perhaps best suited to chat or meta. $\endgroup$ – wolfies Sep 4 '15 at 15:21
  • 1
    $\begingroup$ Thanks very much for your answer! $\endgroup$ – Frank Lu Sep 4 '15 at 19:07

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.