7
$\begingroup$

I am solving a problem from 17th Kolmogorov competition.

Let $X_1,X_2,\dots,X_{50}$ be IID discrete random variables s.t. the probability of $X_1=-1$ equals $\frac 1 9$, the probability of $X_1=- \frac 2 5$ equals $\frac 4 9$, and the probability of $X_1= \frac 1 5$ equals $\frac 4 9$. What is the probability of $X_1+X_2+\cdots +X_{50}=-17$?

Here is the straightforward approach with Mathematica 13.2 on Windows 10.

data = {-1, -2/5, 1/5};weights = {1, 4, 4};
distr=EmpiricalDistribution[weights->data];
Probability[Sum[x[j], {j, 1, 50}] == -17,
Table[x[j] \[Distributed] distr, {j, 1, 50}]]

There are many many variants for Sum[x[j], {j, 1, 50}], so the execution of the latest command takes forever.

How can I calculate this probability?

$\endgroup$
6
  • $\begingroup$ This is a very interesting post. Upvoted all contributors. Thanks. $\endgroup$ Commented Apr 22, 2023 at 7:11
  • $\begingroup$ @PeterMortensen: I deleted your edit in the title and restored the original title. This is at least the second time of such behavior of you. $\endgroup$
    – user64494
    Commented Apr 24, 2023 at 5:30
  • $\begingroup$ Finally located the origin with solution in Russia. The distances are (3/5, 3/5 ) are equal. Transform X -> 1/3 (5(X+1)) and apply Binomi to (1/9,4/9,4/9) -> {0, 1, 2} $\endgroup$
    – Roland F
    Commented Apr 24, 2023 at 11:21
  • 1
    $\begingroup$ Does this answer your question? Convolve discrete random variables efficiently $\endgroup$
    – wolfies
    Commented Apr 24, 2023 at 17:45
  • $\begingroup$ @wolfes: A probability generating function is not implemented in Mathematica as yet. This is noticed in my comment to an answer of JimB. $\endgroup$
    – user64494
    Commented Apr 24, 2023 at 17:56

6 Answers 6

10
$\begingroup$

Another approach that reduces the number of combinations to consider.

Look at the possible combinations of the frequencies (f1, f2, and f3) of a sample size of 50:

Reduce[{f1 (-1) + f2 (-2/5) + f3 (1/5) == -17, 
  f1 + f2 + f3 == 50, 
  0 <= f1 <= 50,
  0 <= f2 <= 50,
  0 <= f3 <= 50}, 
  {f1, f2, f3}]

(* 0 <= f1 <= 45/2 && f2 == 45 - 2 f1 && f3 == 5 (-17 + f1 + (2 f2)/5) *)

Now construct a table of the possible triplets of frequencies such that the sum of the associated values is -17:

ss = Table[{f1, 45 - 2 f1, 5 (-17 + f1 + (2 (45 - 2 f1))/5)}, {f1, 0, 45/2}]
(* {{0, 45, 5}, {1, 43, 6}, {2, 41, 7}, {3, 39, 8}, {4, 37, 9}, 
    {5, 35, 10}, {6, 33, 11}, {7, 31, 12}, {8, 29, 13}, {9, 27, 14}, 
    {10, 25, 15}, {11, 23, 16}, {12, 21, 17}, {13, 19, 18}, {14, 17, 19},
    {15, 15, 20}, {16, 13, 21}, {17, 11, 22}, {18, 9, 23}, {19, 7, 24},
    {20,  5, 25}, {21, 3, 26}, {22, 1, 27}} *)

The desired probability is found using the multinomial formula:

((50!/(#[[1]]! #[[2]]! #[[3]]!)) (1/9)^#[[1]]*(4/9)^(50 - #[[1]]) & /@ ss) // Total
(* 737971498833254262266529255663310449944821760/
   171792506910670443678820376588540424234035840667 *)
$\endgroup$
1
  • $\begingroup$ I prefer that way which basically coincides with the solution from the competition. $\endgroup$
    – user64494
    Commented Apr 21, 2023 at 8:54
10
$\begingroup$

Just to illustrate use of built-in function: MutlinomialDistribution

sol = {a, b, c} /. 
   Solve[{{-1, -2/5, 1/5}, {1, 1, 1}} . {a, b, c} == {-17, 50} && 
     0 <= a <= 50 && 0 <= b <= 50 && 0 <= c <= 50, {a, b, c}, 
    Integers];
pd = MultinomialDistribution[50, {1/9, 4/9, 4/9}];
Total[Probability[{x, y, z} == #, {x, y, z} \[Distributed] pd] & /@ 
  sol]

yields:

737971498833254262266529255663310449944821760/
171792506910670443678820376588540424234035840667

I have voted for @JimB answers: they illustrate how to derive from first principles.

$\endgroup$
1
  • $\begingroup$ +1 The use of MultinomialDistribution with Probability means no messing up when putting in the factorials. $\endgroup$
    – JimB
    Commented Apr 22, 2023 at 4:19
9
$\begingroup$

One way is to use probability generating functions.

pgf = ((1/9) s^(-1) + 4/9 s^(-2/5) + (4/9) s^(1/5))^50
Coefficient[pgf, s^(-17)]
(* 737971498833254262266529255663310449944821760/
171792506910670443678820376588540424234035840667 *)

Edit: The distribution of the sum is found as follows:

distribution = (List @@ Expand[pgf]) /. Times[prob_, Power[s, x_]] -> {x, prob};
ListPlot[distribution, PlotRange -> All]

Distribution of the sum

$\endgroup$
2
  • $\begingroup$ Unfortunately, I don't find a command for pgf in the documentation. $\endgroup$
    – user64494
    Commented Apr 21, 2023 at 8:44
  • $\begingroup$ Neat! and fast. :) $\endgroup$
    – ubpdqn
    Commented Apr 21, 2023 at 11:04
6
$\begingroup$

@ubpdqn 's answer using MultinomialDistribution suggests another approach:

dist = TransformedDistribution[(-1) x1 + (-2/5) x2 + (1/5) x3,
  {x1, x2, x3} \[Distributed] MultinomialDistribution[50, {1/9, 4/9, 4/9}]];
Probability[z == -17, z \[Distributed] dist]]
(* (737971498833254262266529255663310449944821760/
171792506910670443678820376588540424234035840667) *)
$\endgroup$
1
  • 1
    $\begingroup$ glad to see the variety…sometimes TransformedDistribution falter but this is a nice distribution and linear transformation. $\endgroup$
    – ubpdqn
    Commented Apr 23, 2023 at 0:13
6
$\begingroup$

Just to be different (or difficult, nobody really knows), we'll use the CIF.

First transform the generating function so all exponents are integers.

ee = (1/9*t^(-5) + 4/9*t^(-2) + 4/9*t^(1))^(50);

Check that we get the coefficient others have found, using an appropriate transformation of the exponent desired.

In[502]:= Coefficient[ee, t^(-5*17)]

(* Out[502]= \
737971498833254262266529255663310449944821760/\
171792506910670443678820376588540424234035840667 *)

Now prepare it for the CIF. Among other things this involves moving the power in question so that it's -1.

eeb = t^(17*5 - 1)*ee /. t -> Exp[I*theta];

Throw in appropriate factor for the transformation of the differential, and integrate over the unit circle contour to get the residue.

Integrate[Exp[I*theta]*eeb, {theta, 0, 2*Pi}]/(2*Pi)

(* Out[503]= \
737971498833254262266529255663310449944821760/\
171792506910670443678820376588540424234035840667 *)
$\endgroup$
2
  • 1
    $\begingroup$ +1 "different or difficult". I think we do know. ;) $\endgroup$
    – JimB
    Commented Apr 22, 2023 at 17:43
  • $\begingroup$ +1: enjoying the variety on display to crack this humble nut :) $\endgroup$
    – ubpdqn
    Commented Apr 23, 2023 at 0:09
4
$\begingroup$

My solution is as follows.

sol = Solve[-17 == -x - 2/5*y + 1/5*z && {x, y, z} >= 0 && 
x + y + z == 50, {x, y, z}, Integers];
Total[Table[(x + y + z)!/x!/y!/z!*(1/9)^x*(4/9)^(y + z) /. 
sol[[j]], {j, 1, Dimensions[sol][[1]]}]]

737971498833254262266529255663310449944821760/ 171792506910670443678820376588540424234035840667

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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