# ProbabilityDistribution[] function just hangs

I have this discrete probability function that I want to utilize with Mathematica's ProbabilityDistribution[] function

The distribution is:

Pn[0, \[Alpha]_] := 1 - \[Alpha];
Pn[1, \[Alpha]_] := (1 - \[Alpha]) (E^\[Alpha] - 1);
Pn[n_, \[Alpha]_] := (1 - \[Alpha])*
Sum[E^(j*\[Alpha])*(-1)^(n -
j)*(((j*\[Alpha] + n - j)*(j*\[Alpha])^(n - j - 1))/(n - j)!),
{j, 0, n}];

f[k_] := Pn[k - 1, \[Alpha]]/Sum[Pn[j, \[Alpha]], {j, 0, k - 1}] -
Pn[k, \[Alpha]]/Sum[Pn[j, \[Alpha]], {j, 0, k}]


and I would like to use other statistics functions like MomentGeneratingFunction[], ProductDistribution[], etc.

When I type in this:

pmf = ProbabilityDistribution[f[k], {k, 1, 20}];


It just hangs and does nothing...I've waited for over an hour sometimes...I don't even know if I can do this.

Is there a better way to do this? I.e., create a pdf using the function above

• Your f[k] seems quite computationally demanding to me. Have you tried calculating a few values for it on its own, to check whether it is the cause of the hang? Is f[k] already normalized? Commented Feb 19, 2022 at 14:45
• Yes I have - Sum[f[k], {k, 1, 20}] /. {\[Alpha] -> 0.5} computes quickly...
– PiE
Commented Feb 19, 2022 at 14:49
• I am unclear precisely what you mean by normalization, but the pmf does sum to (almost) 1 - even for {k, 1, 20}. Theoretically, the support is $[1,\infty]$
– PiE
Commented Feb 19, 2022 at 14:52
• Evaluating Sum[f[k] /. α -> 7/8, {k, 1, 20}] // Simplify // N[#, 20]& shows that either {k, 1, 20} is not a valid range for the distribution or that the distribution needs to be normalized. Note also that ProbabilityDistribution is not a "pmf" but rather a distribution, call it dist[α]. Then the "pmf" is given by PDF[dist[α], x]. Commented Feb 19, 2022 at 15:41

Clear["Global*"]


Start by preventing Pn from evaluating unless n has a nonnegative integer value.

Pn[0, α_] := 1 - α;
Pn[1, α_] := (1 - α) (E^α - 1);
Pn[n_Integer?NonNegative, α_] := (1 - α)*
Sum[E^(j*α)*(-1)^(n - j)*(((j*α + n - j)*
(j*α)^(n - j - 1))/(n - j)!), {j, 0, n}]

f[k_] := Pn[k - 1, α]/Sum[Pn[j, α], {j, 0, k - 1}] -
Pn[k, α]/Sum[Pn[j, α], {j, 0, k}]


Checking the total probability,

Sum[f[k] /. α -> 7/8, {k, 1, 20}] // Simplify // N[#, 20] &

(* 0.99852321012132571990 *)


Since the total probability for {k, 1, 20} does not equal 1 (and the discrepancy worsens for larger α), the distribution must be normalized.

dist[α_] = ProbabilityDistribution[
f[k], {k, 1, 20, 1},
Assumptions -> {0 < α < 1},
Method -> "Normalize"];


Verifying that the total probability is 1

dist[α] /. ProbabilityDistribution -> Sum // Simplify

(* 1 *)


For symbolic values of α computations are very complicated, e.g.,

LeafCount[Mean[dist[α]] // Simplify] // AbsoluteTiming

(* {1.5569, 3616} *)

LeafCount[MomentGeneratingFunction[dist[α], t] // Simplify] //
AbsoluteTiming

(* {31.216, 5629} *)


For numeric values of α the computations are more responsive

Mean[dist[#]] & /@ Range[0.1, 0.9, 0.1] // AbsoluteTiming

(* {1.21928, {1.10051, 1.20419, 1.31467, 1.43658, 1.57661, 1.74568, 1.96436,
2.27823, 2.77035}} *)

Variance[dist[#]] & /@ Range[0.1, 0.9, 0.1] // AbsoluteTiming

(* {2.97234, {0.101521, 0.212579, 0.345012, 0.51682, 0.760525, 1.14449, 1.83958,
3.37093, 7.15482}} *)

MomentGeneratingFunction[dist[0.5], t] // AbsoluteTiming

(* {0.12198,
E^t (0.606531 + 0.264003 E^t + 0.0910935 E^(2 t) + 0.0274169 E^(3 t) +
0.00784581 E^(4 t) + 0.0022263 E^(5 t) + 0.000632653 E^(6 t) +
0.000180012 E^(7 t) + 0.0000512389 E^(8 t) + 0.0000145858 E^(9 t) +
4.15209*10^-6 E^(10 t) + 1.18197*10^-6 E^(11 t) +
3.36467*10^-7 E^(12 t) + 9.57826*10^-8 E^(13 t) +
2.72672*10^-8 E^(14 t) + 7.75901*10^-9 E^(15 t) +
2.19825*10^-9 E^(16 t) + 6.3983*10^-10 E^(17 t) +
1.56433*10^-10 E^(18 t) - 6.36646*10^-12 E^(19 t))} *)
`
• Thanks so much!
– PiE
Commented Feb 19, 2022 at 17:14