If in a binomial distribution, the Bernoulli trials are independent and have different success probabilities, then it is called Poisson Binomial Distribution. Such a question has been previously answered here and here.

How can I do a similar analysis in the case of a multinomial distribution? For instance, if a $k$-sided die is thrown $n$ times and the probabilities of each side showing up changes every time instead of being fixed (as in the case of regular multinomial distribution), how can I calculate the probability mass function of such a distribution? We assume that we have access to $mathbbp_i_1^n$ where $mathbbp_i$ is a vector of length $k$ denoted the probability of each of the $k$ sides showing up in the $i^th$ trial.

Note: I have asked this question on stats.stackexchange as well, but I feel it is more pertinent here.

Preliminary (TL;DR)

Background

In his 1991 publication, Norman C. Beaulieu answered your question w/ what he dubbed, the generalized multinomial distribution (GMD). My explanation will focus on the GMD's utility.

Notation

• # categories $= c$.


• # trials $= t$.


• Random vector $= X = left[beginarrayccccX_1&X_2&cdots&X_cendarrayright]^T$.


• Category responses after $t$ trials vector $= x = left[beginarrayccccx_1&x_2&cdots&x_cendarrayright]^T$.
  
  
  • 
    $sum_k = 1^c x_k = t$.
  
  


• Probability of category response during trial matrix $= p = left[beginarraycccc p_1,1 & p_1,2 & cdots & p_1,c \
  p_2,1 & p_2,2 & cdots & p_2,c \
  vdots & vdots & ddots & vdots \
  p_t,1 & p_t,2 & cdots & p_t,c
  endarrayright]$.


• Pmf of $X = Pleft[X = xright]$.


• 
  $[c] = left1, 2, cdots, cright$.


• 
  Multiset of $[c] = ([c], m) = left1^m(1), 2^m(2), cdots, c^m(c)right$.
  
  
  • 
    $m(i) = x_i$.
  
  


• 
  Permutations of $([c], m) = mathfrakS_([c], m)$.
  
  
  • 
    $cardleft(mathfrakS_([c], m)right) = left(m(1), m(2), cdots, m(c)right)!$.
  
  

Pmf of GMD

$$Pleft[X = xright] = sum_mathfraks in mathfrakS_([c], m) leftprod_k = 1^t leftp_k,mathfraks_krightright$$

So far, I've identified it as being the superclass of 7 distributions! Namely...

• Bernoulli distribution.


• Uniform distribution.


• Categorical distribution.


• Binomial distribution.


• Multinomial distribution.


• Poisson's binomial distribution.


• Generalized multinomial distribution (if your definition of superclass allows self-inclusion).

Examples

Games

• g1: A 2 sided die is simulated using a fair standard die by assigning faces w/ pips 1 through 3 & 4 through 6 to sides 1 & 2, respectively. The die is biased by etching micro holes into faces w/ pips 1 through 3 s.t. $p_1 = 12/30$ & $p_2 = 18/30$. The 2 sided die is tossed 1 time & the category responses are recorded.


• g2: Same as g1, accept w/ ideal standard die, i.e., $p_1 = p_2 = cdots = p_6 = 5/30$.


• g3: Same as g1, accept w/ standard die, i.e., $p_1 = p_2 = p_3 = 4/30$ & $p_4 = p_5 = p_6 = 6/30$.


• g4: Same as g1, accept die is tossed 7 times.


• g5: Same as g3, accept die is tossed 7 times.


• g6: Same as g4, accept the micro holes are filled w/ $0.07$ kg of a material, which evaporates @ $0.01$ kg/s upon being sprayed w/ an activator, s.t. $p_1 = p_2 = 15/30$ for the 1st toss. Immediately after being sprayed, category responses are recorded every second.


• g7: Same as g6, accept w/ standard die, i.e., $p_1 = p_2 = cdots = p_6 = 5/30$ for the 1st toss.

Questions

• q1: Find pmf & evaluate when $x = left[beginarraycc0&1endarrayright]^T$.


• q2: Find pmf & evaluate when $x = left[beginarraycccccc0&1&0&0&0&0endarrayright]^T$.


• q3: q2.


• q4: Find pmf & evaluate when $x = left[beginarraycc2&5endarrayright]^T$.


• q5: Find pmf & evaluate when $x = left[beginarraycccccc0&2&1&1&0&3endarrayright]^T$.


• q6: q4.


• q7: q5.

Answers w/o knowledge of GMD

• a1: $X$ ~ Bernoulli distribution.
  
  
  • 
    $Pleft[X = xright] = t!prod_k = 1^c fracp_k^kk! = 1!prod_k = 1^2 fracp_k^kk! = frac1!(12/30)^0(18/30)^10!1!$
    $Longrightarrow Pleft[X = xright] = 3/5$.
  
  


• a2: $X$ ~ Uniform distribution.
  
  
  • 
    $Pleft[X = xright] = t!prod_k = 1^c fracp_k^kk! = 1!prod_k = 1^6 fracp_k^kk! = frac1!(5/30)^0 + 1 + 0 + 0 + 0 + 00!1!0!0!0!0!$
    $Longrightarrow Pleft[X = xright] = 1/6$.
  
  


• a3: $X$ ~ Categorical distribution.
  
  
  • 
    $Pleft[X = xright] = t!prod_k = 1^c fracp_k^kk! = 1!prod_k = 1^6 fracp_k^kk! = frac1!(4/30)^0 + 1 + 0(6/30)^0 + 0 + 00!1!0!0!0!0!$
    $Longrightarrow Pleft[X = xright] = 2/15$.
  
  


• a4: $X$ ~ Binomial distribution.
  
  
  • 
    $Pleft[X = xright] = t!prod_k = 1^c fracp_k^kk! = 7!prod_k = 1^2 fracp_k^kk! = frac7!(12/30)^2(18/30)^52!5!$
    $Longrightarrow Pleft[X = xright] = 20412/78125$.
  
  


• a5: $X$ ~ Multinomial distribution.
  
  
  • 
    $Pleft[X = xright] = t!prod_k = 1^c fracp_k^kk! = 7!prod_k = 1^6 fracp_k^kk! =
    frac7!(4/30)^0 + 2 + 1(6/30)^1 + 0 + 30!2!1!1!0!3!$
    $Longrightarrow Pleft[X = xright] = 224/140625$.
  
  


• a6: $X$ ~ Poisson's binomial distribution.
  
  
  • 
    $Pleft[left[beginarrayccX_1&X_2endarrayright]^T = left[beginarrayccx_1&x_2endarrayright]^Tright] = Pleft[X_1 = x_1, X_2 = x_2right] = Pleft[X_1 = x_1right] = Pleft[X_2 = x_2right]$.
  
  
  • 
    $p_1$ & $p_2$ are vectors now: $p_1 = left[beginarrayccccp_1_1&p_1_2&cdots&p_1_tendarrayright]^T, p_2 = left[beginarrayccccp_2_1&p_2_2&cdots&p_2_tendarrayright]^T$.
  
  
  • 
    $Pleft[X_2 = x_2right] = frac1t + 1sum_i = 0^t leftexpleft(frac-j2pi i x_2t + 1right) prod_k = 1^t leftp_2_kleft(expleft(fracj2pi it + 1right) - 1right) + 1rightright$
    $= frac18sum_i = 0^7 leftexpleft(frac-j5pi i4right) prod_k = 1^7 leftleft(frac0.5k + 14.530right)left(expleft(fracjpi i4right) - 1right) + 1rightright$
    $Longrightarrow Pleft[X_2 = 5right] = 308327/1440000$.
  
  


• a7: $X$ ~ Generalized multinomial distribution.
  
  
  • ???
  
  

Answers w/ Knowledge of GMD

• a1: $X$ ~ Bernoulli distribution.
  
  
  • 
    $p = left[beginarraycfrac1230&frac1830endarrayright]$.
  
  
  • 
    $mathfrakS_([2], m) = leftleft(2right)right$.
  
  


• a2: $X$ ~ Uniform distribution.
  
  
  • 
    $p = left[beginarraycfrac530&frac530&frac530&frac530&frac530&frac530endarrayright]$.
  
  
  • 
    $mathfrakS_([6], m) = leftleft(2right)right$.
  
  


• a3: $X$ ~ Categorical distribution.
  
  
  • 
    $p = left[beginarraycfrac430&frac430&frac430&frac630&frac630&frac630endarrayright]$.
  
  
  • 
    $mathfrakS_([6], m) = leftleft(2right)right$.
  
  


• a4: $X$ ~ Binomial distribution.
  
  
  • 
    $p = left[beginarraycc
    frac1230&frac1830 \
    frac1230&frac1830 \
    frac1230&frac1830 \
    frac1230&frac1830 \
    frac1230&frac1830 \
    frac1230&frac1830 \
    frac1230&frac1830
    endarrayright]$.
  
  
  • 
    $mathfrakS_([2], m) = leftleft(1,1,2,2,2,2,2right), ldots, left(2,2,2,2,2,1,1right)right$.
  
  


• a5: $X$ ~ Multinomial distribution.
  
  
  • 
    $p = left[beginarraycccccc
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630 \
    frac430&frac430&frac430&frac630&frac630&frac630
    endarrayright]$.
  
  
  • 
    $mathfrakS_([6], m) = leftleft(2,2,3,4,6,6,6right), ldots, left(6,6,6,4,3,2,2right)right$.
  
  


• a6: $X$ ~ Poisson's binomial distribution.
  
  
  • 
    $p = left[beginarraycc
    frac1530&frac1530 \
    frac14.530&frac15.530 \
    frac1430&frac1630 \
    frac13.530&frac16.530 \
    frac1330&frac1730 \
    frac12.530&frac17.530 \
    frac1230&frac1830
    endarrayright]$.
  
  
  • 
    $mathfrakS_([2], m) = leftleft(1,1,2,2,2,2,2right), ldots, left(2,2,2,2,2,1,1right)right$.
  
  


• a7: $X$ ~ Generalized multinomial distribution.
  
  
  • 
    $p = left[beginarraycccccc
    frac530&frac530&frac530&frac530&frac530&frac530 \
    frac4.8overline330&frac4.8overline330&frac4.8overline330
    &frac5.1overline630&frac5.1overline630&frac5.1overline630 \
    frac4.overline630&frac4.overline630&frac4.overline630
    &frac5.overline330&frac5.overline330&frac5.overline330 \
    frac4.530&frac4.530&frac4.530
    &frac5.530&frac5.530&frac5.530 \
    frac4.overline330&frac4.overline330&frac4.overline330
    &frac5.overline630&frac5.overline630&frac5.overline630 \
    frac4.1overline630&frac4.1overline630&frac4.1overline630
    &frac5.8overline330&frac5.8overline330&frac5.8overline330 \
    frac430&frac430&frac430&frac630&frac630&frac630
    endarrayright]$.
  
  
  • 
    $mathfrakS_([6], m) = leftleft(2,2,3,4,6,6,6right), ldots, left(6,6,6,4,3,2,2right)right$.
  
  
  • 
    $Pleft[X = xright] = 59251/36905625$.
  
  

Final Words

I know my answer was very long (& went far beyond what OP asked for) but this had been flying around inside my head for quite some time & this q seemed like the most suitable landing strip.

I performed the last 6 calculations using the function gmdPmf (which I defined in Mathematica)...

(* GENERALIZED MULTINOMIAL DISTRIBUTION (GMD) *)
(* Note: mXn = # rows X # columns. *)
gmdPmf[
 x_ (* Responses of category j, after t trials have taken place. *),
 p_ (* Matrix (tXm) holds p_trial i, category j = P["Response of trial i is category j"]. *)
] := Module[t, c, ⦋c⦌, allRPs, desiredRPs, count = 0, sum = 0, product = 1,
 t = Total[x]; (* # trials. *)
 c = Length[x]; (* # categories. *)
 ⦋c⦌ = Range[c]; (* Categories. *)
 allRPs = Tuples[⦋c⦌,t]; (* Matrix (c^tXt) holds all the response patterns given that t trials have occurred. *)
 desiredRPs = ; (* Matrix ((x_1,x_2,...,x_c) !Xt) holds the desired response patterns; subset of allRPs wrt n. *)

 For[i = 1, i <= Length[allRPs], i++,
 For[j = 1, j <= c, j++, If[Count[allRPs[[i]],⦋c⦌[[j]]] == x[[j]], count++];];
 If[count == c, AppendTo[desiredRPs, allRPs[[i]]]];
 count = 0;
 ];

 For[i = 1, i <= Length[desiredRPs], i++, 
 For[j = 1, j <= t, j++, product *= (p[[j]][[desiredRPs[[i]][[j]]]]);];
 sum += product;
 product = 1;
 ];

 sum
];

(* ANSWERS *)
Print["a1: P[X = x] = ", gmdPmf[0, 1, 12/30, 18/30], "."];
Print["a2: P[X = x] = ", gmdPmf[0,1, 0, 0, 0, 0, 5/30, 5/30, 5/30, 5/30, 5/30, 5/30], "."];
Print["a3: P[X = x] = ", gmdPmf[0,1, 0, 0, 0, 0, 4/30, 4/30, 4/30, 6/30, 6/30, 6/30], "."];
Print["a4: P[X = x] = ", gmdPmf[2, 5, ArrayFlatten[ConstantArray[12/30, 18/30, 7, 1]]], "."];
Print["a5: P[X = x] = ", gmdPmf[0, 2, 1, 1, 0, 3, ArrayFlatten[ConstantArray[4/30, 4/30, 4/30, 6/30, 6/30, 6/30, 7, 1]]], "."];
p = ; For[i = 1, i <= 7, i++, l = ((31/2) - (1/2)*i)/30; r = ((29/2) + (1/2)*i)/30; AppendTo[p,l,r];]; Print["a6: P[X = x] = ", gmdPmf[2, 5, p], "."];
p = ; For[i = 1, i <= 7, i++, l = ((31/6) - (1/6)*i)/30; r = ((29/6) + (1/6)*i)/30; AppendTo[p,l,l,l,r,r,r];]; Print["a7: P[X = x] = ", gmdPmf[0, 2, 1, 1, 0, 3, p], "."];

Clear[gmdPmf];

Please edit, if you know of any ways to make it shorter/faster. Congrats, if you made it to the end! (:

I have come across very few resources on the topic. Such a problem is called Generalized Poisson's distribution or GPB. I am listing a few of them here to help others.

1. An old paper describing the calculation of its approximate pdf.


2. A 2018 paper describing an algorithm to compute its PDF using Fourier Transforms.


3. An R package implementing the same.