My problem is the following: given a list of probabilities $p_0, p_1, \dots, p_M$ and an integer $n$, I need to compute a function of a real variable $x$

$$W(x) = \sum_{m_1 = 0}^M \sum_{m_2 = 0}^M \dots \sum_{m_n = 0}^M \left( \prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}\right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x), $$ where $L_i (x)$ is the Laguerre polynomial LaguerreL[i,x].

My current (inefficient) solution is this:

w[x_, n_, probs_] :=   (* probs = {p0, p1, p2, .., pM} *)
Module[ {mM = Length@probs -1 , mj },
    { sumIt = Sequence@@Table[ { mj[i], 0 , mM}, {i, 1, n}]},
    Module[ { sumMJ = Sum[ mj[j], {j,1, n}]},
    Product[ probs[[ mj[j]+1]]/ ((mj[j])!), { j , 1, n}]
    sumMJ! ( -1)^sumMJ / n^sumMJ LaguerreL[ sumMJ , x ]], 

It takes quite a long time to compute, e.g. with $n=9$ and $M = 3$ it returns the result after about 2 minutes. I guess this is because the approach is quite procedural. I would like to derive this function for the values $n \leq 40$ and $M \leq 5$, which I am afraid would take a few weeks.

Therefore, my question is, is it possible to compute this function more efficiently in Mathematica?


1 Answer 1


2nd Update: I've added the rationale I used to obtain the more compact sum.

Update: This is now (hopefully) an answer. At end of the case where $M=1$ is the general case.

Original "answer":

This is an extended comment. By setting $M$ to some value and then looking for a pattern for various values of $n$, one might find a quicker formula.

For example, consider $M=1$. Looking at the coefficients for a few values of $n$ and using FindSequenceFunction, I get the following:

c[n_] := DifferenceRoot[
  Function[{y, i}, {(-(n + 1) + i) y[i] + n y[1 + i] == 0,  y[1] == 1}]][#] & /@ Range[n + 1]

w2[x_, n_] := Table[p[0]^(n - i) (-1)^i p[1]^i LaguerreL[i, x], {i, 0, n}] . c[n]

AbsoluteTiming[w2[x, 15];]
(* {0.013407, Null} *)
AbsoluteTiming[w[x, 15, {p[0], p[1]}];]
(* {8.23932, Null} *)

AbsoluteTiming[w2[x, 40];]
(* {0.0418239, Null} *)
AbsoluteTiming[w[x, 40, {p[1], p[2]}];]
(* Still waiting after several minutes *)

General case:

w3[x_, n_, m_] := Module[{t, r},
  t = Flatten[Permutations[#] & /@ (PadRight[#, m + 1] & /@IntegerPartitions[n, m + 1]), 1];
  r = Range[0, m];
  c1 = (n!/Times @@ (#!)) (1/Times @@ (r!^#)) (r . #)!/(-n)^(r . #) & /@ t;
  c2 = Product[p[i]^#[[i + 1]], {i, 0, m}] & /@ t;
  c3 = LaguerreL[r . #, x] & /@ t;
  Total[c1 c2 c3]

AbsoluteTiming[w[x, 3, {p[0], p[1], p[2]}]]

Result with OPs code n=3 and M=2

AbsoluteTiming[w3[x, 3, 2]]

Faster code

Check on if the results are the same:

w[x, 3, {p[0], p[1], p[2]}] - w3[x, 3, 2]
(* 0 *)

A more lengthy calculation:

AbsoluteTiming[w3[x, 40, 4];]
(* {251.448, Null} *)

Rationale for the more compact summation:

We have the following sum:

$$W(x) = \sum_{m_1 = 0}^M \sum_{m_2 = 0}^M \dots \sum_{m_n = 0}^M \left( \prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}\right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)$$

Because we can write

$$\prod_{j=1}^n \frac{ p_{m_j}}{ m_j !}=\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j !$$

where $\sum_{i=0}^M k_i=n$, for each combination of the $m_i$ indices the $n(M+1)$ terms are of the form

$$\left(\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j ! \right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)$$

For example, suppose $M=2$, $n=5$, $m_1=0$, $m_2=1$, $m_3=1$, $m_4=2$, and $m_5=1$. The associated term is

$$\left(\prod_{i=0}^M p_i^{k_i}/\prod_{j=1}^n m_j ! \right) \frac{ \left(\sum_{i=1}^n m_i\right)! }{ ( - n)^{\sum_{i=1}^n m_i} } L_{\sum_{i=1}^n m_i } (x)=p_0 p_1^3 p_2 \frac{5!}{0! (1!)^3 2! (-5)^5}L_5 (x)$$

which simplifies to

$$-\frac{12}{625}p_0 p_1^3 p_2 L_5 (x)$$

There are 19 other combinations of the $m_i$ values that result in the same exact term. How to calculate that multiplier? Of the 5 values of $m_i$ there is one 0, three 1's, and one 2. So 1, 3, and 1 are the powers of $p_0$, $p_1$, and $p_2$, respectively. The desired count is given by


The coefficient associated with $p_0 p_1^3 p_2 L_5 (x)$ is then


So (again, for this example) we just need to generate all of the combinations of $p_0^{k_0} p_1^{k_1} p_2^{k_2}$ where $k_i \geq 0$ and $k_0+k_1+k_2=n$, construct the associated term, and then multiply by the number of combinations of the $m_i$ that will generate that term (which is $\frac{n!}{k_0!k_1!k_2!}$).

  • $\begingroup$ Thank you very much, this is very impressive. I would absolutely love to learn more about the way you were able to come up with this solution. $\endgroup$
    – And R
    Commented Dec 7, 2023 at 8:48

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.