How to find this probability?

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?

• This is a very interesting post. Upvoted all contributors. Thanks. Commented Apr 22, 2023 at 7:11
• @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. Commented Apr 24, 2023 at 5:30
• 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} Commented Apr 24, 2023 at 11:21
• Does this answer your question? Convolve discrete random variables efficiently Commented Apr 24, 2023 at 17:45
• @wolfes: A probability generating function is not implemented in Mathematica as yet. This is noticed in my comment to an answer of JimB. Commented Apr 24, 2023 at 17:56

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 *)

• I prefer that way which basically coincides with the solution from the competition. Commented Apr 21, 2023 at 8:54

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.

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

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]


• Unfortunately, I don't find a command for pgf in the documentation. Commented Apr 21, 2023 at 8:44
• Neat! and fast. :) Commented Apr 21, 2023 at 11:04

@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) *)

• glad to see the variety…sometimes TransformedDistribution falter but this is a nice distribution and linear transformation. Commented Apr 23, 2023 at 0:13

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 *)

• +1 "different or difficult". I think we do know. ;)
– JimB
Commented Apr 22, 2023 at 17:43
• +1: enjoying the variety on display to crack this humble nut :) Commented Apr 23, 2023 at 0:09

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