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?

Improve this question
$\endgroup$
6
  • $\begingroup$ This is a very interesting post. Upvoted all contributors. Thanks. $\endgroup$
    – Romke Bontekoe
    Apr 22 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
    Apr 24 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
    Apr 24 at 11:21
  • 1
    $\begingroup$ Does this answer your question? Convolve discrete random variables efficiently $\endgroup$
    – wolfies
    Apr 24 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
    Apr 24 at 17:56

6 Answers 6

Reset to default
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 *)
Improve this answer
$\endgroup$
1
  • $\begingroup$ I prefer that way which basically coincides with the solution from the competition. $\endgroup$
    – user64494
    Apr 21 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.

Improve this answer
$\endgroup$
1
  • $\begingroup$ +1 The use of MultinomialDistribution with Probability means no messing up when putting in the factorials. $\endgroup$
    – JimB
    Apr 22 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

Improve this answer
$\endgroup$
2
  • $\begingroup$ Unfortunately, I don't find a command for pgf in the documentation. $\endgroup$
    – user64494
    Apr 21 at 8:44
  • $\begingroup$ Neat! and fast. :) $\endgroup$
    – ubpdqn
    Apr 21 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) *)
Improve this answer
$\endgroup$
1
  • 1
    $\begingroup$ glad to see the variety…sometimes TransformedDistribution falter but this is a nice distribution and linear transformation. $\endgroup$
    – ubpdqn
    Apr 23 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 *)
Improve this answer
$\endgroup$
2
  • 1
    $\begingroup$ +1 "different or difficult". I think we do know. ;) $\endgroup$
    – JimB
    Apr 22 at 17:43
  • $\begingroup$ +1: enjoying the variety on display to crack this humble nut :) $\endgroup$
    – ubpdqn
    Apr 23 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

Improve this answer
$\endgroup$

Your Answer

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

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