api.moment

moment_tensortrans(A, T, mode)

B = moment_tensortrans(A, T, mode)

Applies a conversion matrix along one dimension of a joint moment array.

This is the mode product of the array with the matrix: every fibre of A along the given dimension is replaced by T times that fibre. Applying it once per dimension realises the Kronecker product of the univariate conversions, which is the structure of every separable edge of the house of moments.

Input:: A: joint moment array of size (n_1+1)x…x(n_d+1) T: (n_mode+1)x(n_mode+1) conversion matrix mode: dimension to transform (1 <= mode <= d)
Output:: B: array of the same size as A

Example

B = moment_tensortrans(A, moment_stirling1(size(A,1)-1), 1);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_tensorsize(A)

sz = moment_tensorsize(A)

Size of a joint moment array with the trailing singleton dimensions removed.

MATLAB cannot represent a trailing singleton dimension, so a column vector of n+1 elements and an (n+1)x1 joint moment array of a degenerate second class are the same object. The convention taken throughout the joint conversions is the first one: trailing singletons are dropped, so that a column vector is treated as the d = 1 case.

Input:

A: joint moment array

Output:

sz: row vector of the extents, with the trailing singleton dimensions: removed

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_tail_from_binomial(b)

t = moment_tail_from_binomial(b)

Converts the binomial moments of a nonnegative integer random variable into its survival (tail) probabilities,

t_m = sum_{j>=m} (-1)^(j-m) * nchoosek(j-1,m-1) * b_j, m >= 1

with t_0 = 1. Inverse of moment_binomial_from_tail. The inversion is exact on the finite box supplied, the matrix being unit upper triangular, but it reconstructs the true tail only if the binomial moments were themselves those of a law supported on 0,…,n.

Input:: b: vector of length n+1 holding b_0,…,b_n
Output:: t: vector of length n+1 holding t_0,…,t_n, with the same orientation as b

Example

t = moment_tail_from_binomial(moment_binomial_from_tail([1, 1, 1, 0]));

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_raw_from_cumulant(kappa)

m = moment_raw_from_cumulant(kappa)

Converts the cumulants kappa_n of a random variable X into its power (raw) moments m_n = E[X^n], by running the exponential-formula recursion forward,

m_n = sum_{k=1}^{n} nchoosek(n-1,k-1) * kappa_k * m_(n-k)

with m_0 = 1. Equivalently m_n = sum_{pi in P(n)} prod_{B in pi} kappa_|B| over the set partitions of {1,…,n}. Inverse of moment_cumulant_from_raw.

Input:

kappa: vector of length n+1 holding kappa_0,…,kappa_n, i.e. kappa(i) is: the cumulant of order i-1. Element 1 is ignored, since kappa_0 = 0 carries no information

Output:

m: vector of length n+1 holding m_0,…,m_n, with the same orientation as: kappa and m(1) = 1

Example

m = moment_raw_from_cumulant(moment_cumulant_from_raw([1, 2, 6, 22]));

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_jointtrans(A, edge)

B = moment_jointtrans(A, edge)

Applies the conversion matrix of one edge of the house of moments along every dimension of a joint moment array. This is the separable (Kronecker product) form shared by all the joint conversions except the cumulant and the central ones.

Input:: A: joint moment array of size (n_1+1)x…x(n_d+1) edge: edge label accepted by moment_housematrix
Output:: B: array of the same size as A holding the converted joint moments

Example

f = moment_jointtrans(m, ‘factorial_from_raw’);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_upfactorial_from_raw(m)

fp = moment_joint_upfactorial_from_raw(m)

Converts joint power (raw) moments into the joint upward-factorial moments f+_(i_1,…,i_d) = E[prod_j N_j(N_j+1)…(N_j+i_j-1)], by the Stirling cycle numbers along every dimension.

Input:

m: array of size (n_1+1)x…x(n_d+1) holding the joint raw: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

fp: array of the same size holding the joint upward-factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_upfactorial_from_negbinomial(bm)

fp = moment_joint_upfactorial_from_negbinomial(bm)

Converts joint negative-binomial moments into joint upward-factorial moments. Inverse of moment_joint_negbinomial_from_upfactorial.

Input:

bm: array of size (n_1+1)x…x(n_d+1) holding the joint negative-binomial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

fp: array of the same size holding the joint upward-factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_upfactorial_from_factorial(f)

fp = moment_joint_upfactorial_from_factorial(f)

Converts joint factorial moments into joint upward-factorial moments, by the Lah numbers along every dimension.

Input:

f: array of size (n_1+1)x…x(n_d+1) holding the joint factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

fp: array of the same size holding the joint upward-factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_tail_from_binomial(b)

t = moment_joint_tail_from_binomial(b)

Converts the joint binomial moments of a nonnegative integer random vector into its joint survival probabilities. Inverse of moment_joint_binomial_from_tail.

Input:: b: array of size (n_1+1)x…x(n_d+1) holding the joint binomial moments
Output:: t: array of the same size holding the joint survival probabilities

Example

t = moment_joint_tail_from_binomial(moment_joint_binomial_from_tail(t0));

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_raw_from_upfactorial(fp)

m = moment_joint_raw_from_upfactorial(fp)

Converts joint upward-factorial moments into joint power (raw) moments, by the signed Stirling numbers of the second kind along every dimension. Inverse of moment_joint_upfactorial_from_raw.

Input:

fp: array of size (n_1+1)x…x(n_d+1) holding the joint upward-factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

m: array of the same size holding the joint raw moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_raw_from_factorial(f)

m = moment_joint_raw_from_factorial(f)

Converts joint factorial moments into joint power (raw) moments, by the Stirling numbers of the second kind along every dimension. Inverse of moment_joint_factorial_from_raw.

Input:

f: array of size (n_1+1)x…x(n_d+1) holding the joint factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

m: array of the same size holding the joint raw moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_raw_from_cumulant(kappa)

m = moment_joint_raw_from_cumulant(kappa)

Converts the joint cumulants of a random vector into its joint power (raw) moments, by running the multivariate exponential-formula recursion forward,

m_a = sum_(0<b<=a) prod_l nchoosek(a_l-[l=j], b_l-[l=j]) kappa_b m_(a-b)

with m_0 = 1 and j the first dimension in which a is nonzero. Inverse of moment_joint_cumulant_from_raw.

Input:

kappa: array of size (n_1+1)x…x(n_d+1) holding the joint cumulants;: element 1 is ignored

Output:

m: array of the same size holding the joint power moments, element 1 being: 1

Example

m = moment_joint_raw_from_cumulant(moment_joint_cumulant_from_raw(m0));

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_joint_raw_from_central(mc, mu)

m = moment_joint_raw_from_central(mc, mu)

Converts the joint central moments of a random vector into the joint power (raw) moments, by the multi-index binomial theorem,

m_(i) = sum_(k<=i) prod_j nchoosek(i_j,k_j) mu_j^(i_j-k_j) * mc_(k)

The mean vector must be supplied separately, since the first-order central moments are zero and carry no information on it. Inverse of moment_joint_central_from_raw.

Input:: mc: array of size (n_1+1)x…x(n_d+1) holding the joint central moments mu: vector of length d holding the means E[N_1],…,E[N_d]
Output:: m: array of the same size as mc holding the joint power moments

Example

m = moment_joint_raw_from_central(moment_joint_central_from_raw(m0), mu);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

moment_joint_negbinomial_from_upfactorial(fp)

bm = moment_joint_negbinomial_from_upfactorial(fp)

Converts joint upward-factorial moments into the joint negative-binomial moments b-_(i_1,…,i_d) = E[prod_j nchoosek(N_j+i_j-1,i_j)] = f+_(i) / prod_j (i_j!).

Input:

fp: array of size (n_1+1)x…x(n_d+1) holding the joint upward-factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

bm: array of the same size holding the joint negative-binomial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_negbinomial_from_binomial(b)

bm = moment_joint_negbinomial_from_binomial(b)

Converts joint binomial moments into joint negative-binomial moments, by the shifted binomial transform nchoosek(i-1,k-1) along every dimension.

Input:

b: array of size (n_1+1)x…x(n_d+1) holding the joint binomial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

bm: array of the same size holding the joint negative-binomial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_marking(f, p, dims)

F = moment_joint_marking(f, p, dims)

Joint factorial moments of the per-class counts under multinomial marking.

If a count N is marked independently, every event receiving class j with probability p_j, then the per-class counts (N_1,…,N_d) have joint factorial moments

E[prod_j (N_j)_(a_j)] = (prod_j p_j^(a_j)) * f_(|a|)

where f is the factorial moment sequence of the aggregate count N and |a| = a_1+…+a_d. This is the counting-process counterpart of the marking (class-splitting) formulas of the M3A fitters, and it is exact for the per-class counts of a MAP marked in this i.i.d. way, in particular for an MMAP whose marking probabilities do not depend on the phase.

Input:

f: vector of length n+1 holding the factorial moments f_0,…,f_n of the: aggregate count

p: vector of length d holding the marking probabilities dims: vector of length d holding the maximum order per class. Their sum

must not exceed n, since an entry of multi-order a consumes the aggregate moment of order |a|

Output:

F: array of size (dims(1)+1)x…x(dims(d)+1) holding the joint factorial: moments of the per-class counts

Example

F = moment_joint_marking([1, 2, 4, 8], [0.3, 0.7], [1, 1]);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_factorial_from_upfactorial(fp)

f = moment_joint_factorial_from_upfactorial(fp)

Converts joint upward-factorial moments into joint factorial moments, by the signed Lah numbers along every dimension.

Input:

fp: array of size (n_1+1)x…x(n_d+1) holding the joint upward-factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

f: array of the same size holding the joint factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_factorial_from_raw(m)

f = moment_joint_factorial_from_raw(m)

Converts the joint power (raw) moments m_(i_1,…,i_d) = E[prod_j N_j^(i_j)] of a random vector (N_1,…,N_d) into the joint factorial moments f_(i_1,…,i_d) = E[prod_j (N_j)_(i_j)], where (N)_i = N(N-1)…(N-i+1), by applying the signed Stirling numbers of the first kind separately along every dimension,

f_(i) = sum_(k) prod_j s(i_j,k_j) * m_(k)

The joint conversion is the Kronecker product of the univariate ones, which is what makes the mode-by-mode evaluation legitimate. Only the cumulant and the central conversions are not of this separable form.

Input:

m: array of size (n_1+1)x…x(n_d+1) holding the joint raw: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

f: array of the same size holding the joint factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_factorial_from_factcumulant(kappa)

f = moment_joint_factorial_from_factcumulant(kappa)

Converts the joint factorial cumulants of a discrete random vector into its joint factorial moments. Inverse of moment_joint_factcumulant_from_factorial.

Input:

kappa: array of size (n_1+1)x…x(n_d+1) holding the joint factorial: cumulants; element 1 is ignored

Output:

f: array of the same size holding the joint factorial moments, element 1: being 1

Example

f = moment_joint_factorial_from_factcumulant(kappa);

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_joint_factorial_from_binomial(b)

f = moment_joint_factorial_from_binomial(b)

Converts joint binomial moments into joint factorial moments. Inverse of moment_joint_binomial_from_factorial.

Input:

b: array of size (n_1+1)x…x(n_d+1) holding the joint binomial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

f: array of the same size holding the joint factorial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_factcumulant_from_factorial(f)

kappa = moment_joint_factcumulant_from_factorial(f)

Converts the joint factorial moments of a discrete random vector into its joint factorial cumulants, the coefficients of the logarithm of the joint probability generating function expanded about z = (1,…,1),

log E[prod_j z_j^(N_j)] = sum_(a ~= 0) kappa_a prod_j (z_j-1)^(a_j) / a_j!

They stand to the joint factorial moments exactly as the joint cumulants stand to the joint power moments, so the same recursion applies. For a multivariate Poisson vector with independent components every joint factorial cumulant of order two or more vanishes; for the per-class counts of a marked MAP they measure the departure from independent Poisson marking.

Input:

f: array of size (n_1+1)x…x(n_d+1) holding the joint factorial moments,: with element 1 equal to 1

Output:

kappa: array of the same size holding the joint factorial cumulants,: element 1 being 0

Example

kappa = moment_joint_factcumulant_from_factorial(f);

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_joint_cumulant_from_raw(m)

kappa = moment_joint_cumulant_from_raw(m)

Converts the joint power (raw) moments of a random vector (N_1,…,N_d) into its joint cumulants, the coefficients of the joint cumulant generating function

log E[exp(s_1 N_1 + … + s_d N_d)] = sum_(a ~= 0) kappa_a prod_j s_j^(a_j) / a_j!

They obey the multivariate exponential formula, equivalently the Leonov-Shiryaev partition formula. With j the first dimension in which the multi-index a is nonzero,

m_a = sum_(0<b<=a) prod_l nchoosek(a_l-[l=j], b_l-[l=j]) kappa_b m_(a-b)

which isolates kappa_a because the b = a term has unit coefficient and m_0 = 1. Unlike every other conversion in the house, this one does not factor into a product of univariate transforms: the cumulant of multi-order (1,1) is the covariance, which mixes the dimensions.

Input:

m: array of size (n_1+1)x…x(n_d+1) holding the joint power moments, with: element 1 equal to 1

Output:

kappa: array of the same size holding the joint cumulants, element 1 being: kappa_0 = 0

Example

kappa = moment_joint_cumulant_from_raw(m); % kappa(2,2) is the covariance

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_joint_central_from_tail(t)

mc = moment_joint_central_from_tail(t)

Joint central moments of a nonnegative integer random vector from its joint survival array, composing the four edges that separate the two vertices: tail -> binomial -> factorial -> raw -> central, with the means read off the raw array.

This is the whole path from a solver that produces survival probabilities (a closed queueing network through its normalizing constants, a CTMC through its stationary distribution, a simulator through a histogram) to the covariances and the higher central moments.

Input:

t: array of size (n_1+1)x…x(n_d+1) holding the joint survival: probabilities and covering the support, element 1 being 1

Output:

mc: array of the same size holding the joint central moments. The entry of: multi-order e_j+e_l is the covariance of N_j and N_l

Example

mc = moment_joint_central_from_tail(t);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

moment_joint_central_from_raw_mean(m, mu)

mc = moment_joint_central_from_raw_mean(m, mu)

Converts the joint power (raw) moments of a random vector into the joint central moments about a given mean vector. Same conversion as moment_joint_central_from_raw, with the means supplied rather than read off the array, so that it also applies when the array does not carry the first-order entries.

Input:: m: array of size (n_1+1)x…x(n_d+1) holding the joint power moments mu: vector of length d holding the means E[N_1],…,E[N_d]
Output:: mc: array of the same size as m holding the joint central moments

Example

mc = moment_joint_central_from_raw_mean(m, [1.5, 2.5]);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

moment_joint_central_from_raw(m)

mc = moment_joint_central_from_raw(m)

Converts the joint power (raw) moments of a random vector (N_1,…,N_d) into the joint central moments mc_(i_1,…,i_d) = E[prod_j (N_j - E N_j)^(i_j)].

The conversion is the multi-index binomial theorem, which is again separable but with a different shift per dimension,

mc_(i) = sum_(k<=i) prod_j (-1)^(i_j-k_j) nchoosek(i_j,k_j) mu_j^(i_j-k_j)

m_(k)

The means mu_j = m_(e_j) are read off the array itself, so every dimension must carry at least the first order. The entry of multi-order e_j+e_l is the covariance of N_j and N_l.

Input:

m: array of size (n_1+1)x…x(n_d+1) holding the joint power moments, with: every n_j >= 1

Output:

mc: array of the same size holding the joint central moments

Example

mc = moment_joint_central_from_raw(m);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

moment_joint_binomial_from_tail(t)

b = moment_joint_binomial_from_tail(t)

Converts the joint survival probabilities t_(m_1,…,m_d) = P(N_1 >= m_1, …, N_d >= m_d) of a nonnegative integer random vector into its joint binomial moments,

b_(k) = E[prod_j nchoosek(N_j,k_j)]
= sum_(m>=k) prod_j nchoosek(m_j-1,k_j-1) * t_(m)

The transform is the tensor product of the univariate one, which is what makes the mode-by-mode evaluation legitimate: an entry with k_j = 0 selects m_j = 0, and t_(0,m_2,…) is by construction the marginal survival array of the remaining coordinates. As in the univariate case it is upper triangular, so the array must cover the joint support to be exact; truncating gives lower bounds.

Input:

t: array of size (n_1+1)x…x(n_d+1) holding the joint survival: probabilities, element 1 being 1

Output:

b: array of the same size holding the joint binomial moments

Example

b = moment_joint_binomial_from_tail(t);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_binomial_from_negbinomial(bm)

b = moment_joint_binomial_from_negbinomial(bm)

Converts joint negative-binomial moments into joint binomial moments, by the signed shifted binomial transform along every dimension. Inverse of moment_joint_negbinomial_from_binomial.

Input:

bm: array of size (n_1+1)x…x(n_d+1) holding the joint negative-binomial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

b: array of the same size holding the joint binomial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_binomial_from_factorial(f)

b = moment_joint_binomial_from_factorial(f)

Converts joint factorial moments into the joint binomial moments b_(i_1,…,i_d) = E[prod_j nchoosek(N_j,i_j)] = f_(i) / prod_j (i_j!).

Input:

f: array of size (n_1+1)x…x(n_d+1) holding the joint factorial: moments, element (i_1+1,…,i_d+1) being the moment of multi-order (i_1,…,i_d) and element 1 being 1

Output:

b: array of the same size holding the joint binomial moments

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_joint_aggregate(F)

f = moment_joint_aggregate(F)

Factorial moments of a total count from the joint factorial moments of its parts. For N = N_1+…+N_d the Vandermonde convolution of falling factorials gives

f_n = sum_(|a|=n) (n! / prod_j a_j!) * F_a

which holds for ANY joint law of the parts, marked or not, and is the inverse direction of moment_joint_marking whenever the marking is multinomial. The order reached is limited by the smallest per-class order in F, since the term a = n*e_j must be available for every j.

Input:

F: array of size (n_1+1)x…x(n_d+1) holding the joint factorial moments of: the parts

Output:

f: column vector of length min_j(n_j)+1 holding f_0,…,f_min_j(n_j), the: factorial moments of the total

Example

f = moment_joint_aggregate(moment_joint_marking([1, 2, 4], [0.3, 0.7], [1, 1]));

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_housematrix(edge, n)

T = moment_housematrix(edge, n)

Conversion matrix of one edge of the house of moments.

The edge is returned as a linear map on the moment subspace {m_0 = 1}. Four edges (the Lah pair and the shifted-binomial pair) pin their zeroth output to 1 rather than propagating element 0, so as maps of the whole space they are affine. Here the offset is folded into column 0, which is empty for those edges, making every edge a genuine matrix. On a moment vector, whose element 0 is 1 by definition, the two agree. This is also what makes those edges usable dimension by dimension in the joint conversions.

Input:

edge: one of ‘factorial_from_raw’, ‘raw_from_factorial’,: ‘upfactorial_from_raw’, ‘raw_from_upfactorial’, ‘binomial_from_factorial’, ‘factorial_from_binomial’, ‘negbinomial_from_upfactorial’, ‘upfactorial_from_negbinomial’, ‘factorial_from_upfactorial’, ‘upfactorial_from_factorial’, ‘negbinomial_from_binomial’, ‘binomial_from_negbinomial’, ‘binomial_from_tail’, ‘tail_from_binomial’

n: maximum order of the mode (n >= 0)

Output:

T: (n+1)x(n+1) conversion matrix

Example

T = moment_housematrix(‘factorial_from_raw’, 4);

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_factorial_from_factcumulant(kappa)

f = moment_factorial_from_factcumulant(kappa)

Converts the factorial cumulants of a discrete random variable N into its factorial moments f_n = E[N(N-1)…(N-n+1)], by running the recursion

f_n = sum_{k=1}^{n} nchoosek(n-1,k-1) * kappa_k * f_(n-k)

forward with f_0 = 1. Inverse of moment_factcumulant_from_factorial.

Input:

kappa: vector of length n+1 holding the factorial cumulants of order: 0,…,n. Element 1 is ignored

Output:

f: vector of length n+1 holding f_0,…,f_n, with the same orientation as: kappa and f(1) = 1

Example

f = moment_factorial_from_factcumulant([0, 2, 0, 0]);

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_factcumulant_from_factorial(f)

kappa = moment_factcumulant_from_factorial(f)

Converts the factorial moments f_n = E[N(N-1)…(N-n+1)] of a discrete random variable N into its factorial cumulants, the coefficients of the logarithm of the probability generating function expanded about z = 1,

log E[z^N] = sum_{n>=1} kappa_n (z-1)^n / n!

The factorial cumulants stand to the factorial moments exactly as the cumulants stand to the power moments, so the same recursion applies,

f_n = sum_{k=1}^{n} nchoosek(n-1,k-1) * kappa_k * f_(n-k)

For a Poisson variable of rate lambda all factorial cumulants beyond the first vanish, which makes them the natural measure of departure from Poisson behaviour in the counting process of a MAP.

Input:

f: vector of length n+1 holding f_0,…,f_n, i.e. f(i) is the moment of: order i-1 and f(1) = 1

Output:

kappa: vector of length n+1 holding the factorial cumulants of order: 0,…,n, with the same orientation as f and element 1 equal to 0

Example

kappa = moment_factcumulant_from_factorial([1, 2, 4, 8]);

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_cumulant_from_raw(m)

kappa = moment_cumulant_from_raw(m)

Converts the power (raw) moments m_n = E[X^n] of a random variable X into its cumulants kappa_n, the coefficients of the cumulant generating function log E[exp(sX)] = sum_{n>=1} kappa_n s^n / n!.

The conversion inverts the exponential-formula recursion

m_n = sum_{k=1}^{n} nchoosek(n-1,k-1) * kappa_k * m_(n-k)

equivalently kappa_n = sum_{pi in P(n)} (|pi|-1)! (-1)^(|pi|-1) prod_{B in pi} m_|B| over the set partitions of {1,…,n}. The first cumulants are kappa_1 = m_1, kappa_2 = m_2 - m_1^2 (the variance) and kappa_3 = m_3 - 3 m_1 m_2 + 2 m_1^3 (the third central moment). The conversion is not restricted to discrete random variables.

Input:

m: vector of length n+1 holding m_0,…,m_n, i.e. m(i) is the moment of: order i-1 and m(1) = 1

Output:

kappa: vector of length n+1 holding kappa_0,…,kappa_n, with the same: orientation as m. Element 1 is kappa_0 = 0, the value of the cumulant generating function at the origin, and not m_0 = 1

Example

kappa = moment_cumulant_from_raw([1, 2, 6, 22]);

Reference: V. P. Leonov and A. N. Shiryaev. On a method of calculation of semi-invariants. Theory of Probability and its Applications, 4(3):319-329, 1959.

moment_binotrans(x)

y = moment_binotrans(x)

Binomial transform of the sequence x_0,x_1,…,x_n into y_0,y_1,…,y_n,

y_n = sum_{k=0}^{n} (-1)^(n-k) * nchoosek(n,k) * x_k

Applied to a moment sequence m_i = E[X^i] it returns the moments of the unit downshift, y_i = E[(X-1)^i]. It is not an involution: its inverse is moment_binotransinv, the unsigned transform.

Input:

x: vector of length n+1 holding x_0,…,x_n, i.e. x(i) is the element of: order i-1

Output:

y: vector of length n+1 holding y_0,…,y_n, with the same orientation: as x

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (8).

Example: y = moment_binotrans([1,2,5,15])

moment_binomial_from_tail(t)

b = moment_binomial_from_tail(t)

Converts the survival (tail) probabilities t_m = P(N >= m) of a nonnegative integer random variable into its binomial moments,

b_j = E[nchoosek(N,j)] = sum_{m>=j} nchoosek(m-1,j-1) * t_m, j >= 1

with b_0 = t_0 = 1. Unlike every other edge of the house of moments this transform is UPPER triangular, so it consumes the whole tail: the result is exact only if the sequence covers the support, i.e. t_m = 0 beyond the last element supplied. This is the natural entry point for a closed queueing network, whose queue lengths are bounded by the population and whose joint survival probabilities are ratios of normalizing constants.

Truncating the tail early yields a strict LOWER bound on every b_j, since all the coefficients and all the tail values are nonnegative. The bound is not inherited by the central moments downstream, whose conversion alternates in sign.

Input:

t: vector of length n+1 holding t_0,…,t_n, i.e. t(i) is P(N >= i-1) and: t(1) = 1

Output:

b: vector of length n+1 holding b_0,…,b_n, with the same orientation as t

Example

b = moment_binomial_from_tail([1, 1, 1, 1, 0]); % N = 3 with probability 1

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003.

moment_upfactorial_from_raw(m)

fp = moment_upfactorial_from_raw(m)

Converts the power (raw) moments m_n = E[N^n] of a discrete random variable N into the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] by means of the Stirling cycle numbers,

f_n^+ = sum_{k=0}^{n} sigma(n,k) * m_k

Upward-factorial moments are of use in moment-matching techniques for matrix-geometric and discrete phase-type distributions.

Input:

m: vector of length n+1 holding m_0,…,m_n, i.e. m(i) is the moment of: order i-1 and m(1) = m_0 = 1

Output:

fp: vector of length n+1 holding f_0^+,…,f_n^+, with the same: orientation as m

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: fp = moment_upfactorial_from_raw([1,2,6,22])

moment_upfactorial_from_negbinomial(bm)

fp = moment_upfactorial_from_negbinomial(bm)

Converts the negative-binomial moments b_n^- = E[nchoosek(N+n-1,n)] of a discrete random variable N into the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] via the one-to-one correspondence

f_n^+ = n! * b_n^-

Input:

bm: vector of length n+1 holding b_0^-,…,b_n^-, i.e. bm(i) is the moment: of order i-1 and bm(1) = b_0^- = 1

Output:

fp: vector of length n+1 holding f_0^+,…,f_n^+, with the same: orientation as bm

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (7).

Example: fp = moment_upfactorial_from_negbinomial([1,2,4,22/3])

moment_upfactorial_from_factorial(f)

fp = moment_upfactorial_from_factorial(f)

Converts the factorial moments f_n = E[N(N-1)…(N-n+1)] of a discrete random variable N into the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] by means of the Lah numbers,

f_n^+ = sum_{k=1}^{n} L(n,k) * f_k for n >= 1 f_0^+ = 1

Input:

f: vector of length n+1 holding f_0,…,f_n, i.e. f(i) is the moment of: order i-1 and f(1) = f_0 = 1

Output:

fp: vector of length n+1 holding f_0^+,…,f_n^+, with the same: orientation as f

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: fp = moment_upfactorial_from_factorial([1,2,4,8])

moment_stirlingcycle(n)

sigma = moment_stirlingcycle(n)

Triangle of the Stirling cycle numbers (unsigned Stirling numbers of the first kind), sigma(i,j) = (-1)^(i-j) * s(i,j), obtained from the recursion

sigma(i,j) = (i-1)*sigma(i-1,j) + sigma(i-1,j-1) for j > 0 sigma(0,0) = 1, sigma(i,0) = 0 for i > 0

These numbers are the coefficients that convert power moments into upward-factorial moments.

Input:

n: maximum order (n >= 0)

Output:

sigma: (n+1)x(n+1) matrix with sigma(i+1,j+1) = sigma(i,j) in the 0-based: notation of the reference. Entries with j > i are zero.

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (12).

Example: sigma = moment_stirlingcycle(3)

moment_stirling2(n)

S = moment_stirling2(n)

Triangle of the Stirling numbers of the second kind S(i,j), implicitly defined by the expansion of a power into falling factorials

x^i = sum_{j=0}^{i} S(i,j) x(x-1)(x-2)…(x-j+1)

and computed from the recursion

S(i,j) = j*S(i-1,j) + S(i-1,j-1) for j > 0 S(0,0) = 1, S(i,0) = 0 for i > 0

These numbers are the coefficients that convert factorial moments back into power moments.

Input:

n: maximum order (n >= 0)

Output:

S: (n+1)x(n+1) matrix with S(i+1,j+1) = S(i,j) in the 0-based notation of: the reference. Entries with j > i are zero.

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (11).

Example: S = moment_stirling2(3)

moment_stirling1(n)

s = moment_stirling1(n)

Triangle of the signed Stirling numbers of the first kind s(i,j), defined as the coefficients of x^j in the falling factorial

sum_{j=0}^{i} s(i,j) x^j = x(x-1)(x-2)…(x-i+1)

These numbers are the coefficients that convert power moments into factorial moments. They relate to the Stirling cycle numbers via s(i,j) = (-1)^(i-j) * sigma(i,j).

Input:

n: maximum order (n >= 0)

Output:

s: (n+1)x(n+1) matrix with s(i+1,j+1) = s(i,j) in the 0-based notation of: the reference. Entries with j > i are zero.

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (10) and eq. (12).

Example: s = moment_stirling1(3)

moment_raw_from_upfactorial(fp)

m = moment_raw_from_upfactorial(fp)

Converts the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] of a discrete random variable N into the power (raw) moments m_n = E[N^n] by means of the signed Stirling numbers of the second kind,

m_n = sum_{k=0}^{n} (-1)^(n-k) * S(n,k) * f_k^+

Input:

fp: vector of length n+1 holding f_0^+,…,f_n^+, i.e. fp(i) is the moment: of order i-1 and fp(1) = f_0^+ = 1

Output:

m: vector of length n+1 holding m_0,…,m_n, with the same orientation: as fp

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: m = moment_raw_from_upfactorial(moment_upfactorial_from_raw([1,2,6,22]))

moment_raw_from_factorial(f)

m = moment_raw_from_factorial(f)

Converts the factorial moments f_n = E[N(N-1)…(N-n+1)] of a discrete random variable N into the power (raw) moments m_n = E[N^n] by means of the Stirling numbers of the second kind,

m_n = sum_{k=0}^{n} S(n,k) * f_k

Input:

f: vector of length n+1 holding f_0,…,f_n, i.e. f(i) is the moment of: order i-1 and f(1) = f_0 = 1

Output:

m: vector of length n+1 holding m_0,…,m_n, with the same orientation: as f

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (13).

Example: m = moment_raw_from_factorial(moment_factorial_from_raw([1,2,6,22]))

moment_raw_from_central(mc, m1)

m = moment_raw_from_central(mc, m1)

Converts the central moments m_n^c = E[(N-m_1)^n] of a random variable N into the power (raw) moments m_n = E[N^n] by means of the inverse binomial transform in the variation that involves the mean m_1,

m_n = sum_{k=0}^{n} nchoosek(n,k) * m_k^c * m_1^(n-k)

The mean must be supplied separately since m_1^c = 0 carries no information on it. The conversion also holds for continuous random variables.

Input:

mc: vector of length n+1 holding m_0^c,…,m_n^c, i.e. mc(i) is the moment: of order i-1 and mc(1) = m_0^c = 1

m1: mean of N

Output:

m: vector of length n+1 holding m_0,…,m_n, with the same orientation: as mc

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: m = moment_raw_from_central(moment_central_from_raw([1,2,6,22]), 2)

moment_negbinomial_from_upfactorial(fp)

bm = moment_negbinomial_from_upfactorial(fp)

Converts the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] of a discrete random variable N into the negative-binomial moments b_n^- = E[nchoosek(N+n-1,n)] via the one-to-one correspondence

b_n^- = f_n^+ / n!

Input:

fp: vector of length n+1 holding f_0^+,…,f_n^+, i.e. fp(i) is the moment: of order i-1 and fp(1) = f_0^+ = 1

Output:

bm: vector of length n+1 holding b_0^-,…,b_n^-, with the same: orientation as fp

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (7).

Example: bm = moment_negbinomial_from_upfactorial([1,2,8,44])

moment_negbinomial_from_binomial(b)

bm = moment_negbinomial_from_binomial(b)

Converts the binomial moments b_n = E[nchoosek(N,n)] of a discrete random variable N into the negative-binomial moments b_n^- = E[nchoosek(N+n-1,n)] by means of the shifted binomial transform

b_n^- = sum_{k=1}^{n} nchoosek(n-1,k-1) * b_k for n >= 1 b_0^- = 1

Input:

b: vector of length n+1 holding b_0,…,b_n, i.e. b(i) is the moment of: order i-1 and b(1) = b_0 = 1

Output:

bm: vector of length n+1 holding b_0^-,…,b_n^-, with the same: orientation as b

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (14).

Example: bm = moment_negbinomial_from_binomial([1,2,2,4/3])

moment_lah(n)

L = moment_lah(n)

Triangle of the Lah numbers L(i,j) = (i!/j!)*nchoosek(i-1,j-1), which link the factorial moments to the upward-factorial moments. The triangle is built from the equivalent recursion

L(i,j) = L(i-1,j-1) + (i+j-1)*L(i-1,j) for j > 0 L(0,0) = 1, L(i,0) = 0 for i > 0

which avoids the overflow of the explicit factorial form for large orders.

Input:

n: maximum order (n >= 0)

Output:

L: (n+1)x(n+1) matrix with L(i+1,j+1) = L(i,j) in the 0-based notation of: the reference. Entries with j > i are zero.

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4. I. Lah. Eine neue Art von Zahlen, ihre Eigenschaften und Anwendung in der mathematischen Statistik. Mitteilungsbl. Math. Statist., 7:203-212, 1955.

Example: L = moment_lah(3)

moment_factorial_from_upfactorial(fp)

f = moment_factorial_from_upfactorial(fp)

Converts the upward-factorial moments f_n^+ = E[N(N+1)…(N+n-1)] of a discrete random variable N into the factorial moments f_n = E[N(N-1)…(N-n+1)] by means of the Lah numbers,

f_n = sum_{k=1}^{n} (-1)^(n-k) * L(n,k) * f_k^+ for n >= 1 f_0 = 1

Input:

fp: vector of length n+1 holding f_0^+,…,f_n^+, i.e. fp(i) is the moment: of order i-1 and fp(1) = f_0^+ = 1

Output:

f: vector of length n+1 holding f_0,…,f_n, with the same orientation: as fp

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: f = moment_factorial_from_upfactorial([1,2,8,44])

moment_factorial_from_raw(m)

f = moment_factorial_from_raw(m)

Converts the power (raw) moments m_n = E[N^n] of a discrete random variable N into the factorial moments f_n = E[N(N-1)…(N-n+1)] by means of the signed Stirling numbers of the first kind,

f_n = sum_{k=0}^{n} s(n,k) * m_k

Input:

m: vector of length n+1 holding m_0,…,m_n, i.e. m(i) is the moment of: order i-1 and m(1) = m_0 = 1

Output:

f: vector of length n+1 holding f_0,…,f_n, with the same orientation: as m

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (13).

Example: f = moment_factorial_from_raw([1,2,6,22])

moment_factorial_from_binomial(b)

f = moment_factorial_from_binomial(b)

Converts the binomial moments b_n = E[nchoosek(N,n)] of a discrete random variable N into the factorial moments f_n = E[N(N-1)…(N-n+1)] via the one-to-one correspondence

f_n = n! * b_n

Input:

b: vector of length n+1 holding b_0,…,b_n, i.e. b(i) is the moment of: order i-1 and b(1) = b_0 = 1

Output:

f: vector of length n+1 holding f_0,…,f_n, with the same orientation: as b

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: f = moment_factorial_from_binomial([1,2,2,4/3])

moment_central_from_raw(m)

mc = moment_central_from_raw(m)

Converts the power (raw) moments m_n = E[N^n] of a random variable N into the central moments m_n^c = E[(N-m_1)^n] by means of the binomial transform in the variation that involves the mean m_1,

m_n^c = sum_{k=0}^{n} (-1)^(n-k) * nchoosek(n,k) * m_k * m_1^(n-k)

The conversion also holds for continuous random variables.

Input:

m: vector of length n+1 holding m_0,…,m_n, i.e. m(i) is the moment of: order i-1 and m(1) = m_0 = 1. At least the mean m_1 must be given, hence n >= 1

Output:

mc: vector of length n+1 holding m_0^c,…,m_n^c, with the same: orientation as m. By construction m_0^c = 1 and m_1^c = 0

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: mc = moment_central_from_raw([1,2,6,22])

moment_binotransinv(y)

x = moment_binotransinv(y)

Inverse binomial transform of the sequence y_0,y_1,…,y_n,

x_n = sum_{k=0}^{n} nchoosek(n,k) * y_k

Input:

y: vector of length n+1 holding y_0,…,y_n, i.e. y(i) is the element of: order i-1

Output:

x: vector of length n+1 holding x_0,…,x_n, with the same orientation: as y

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (9).

Example: x = moment_binotransinv(moment_binotrans([1,2,5,15]))

moment_binomial_from_negbinomial(bm)

b = moment_binomial_from_negbinomial(bm)

Converts the negative-binomial moments b_n^- = E[nchoosek(N+n-1,n)] of a discrete random variable N into the binomial moments b_n = E[nchoosek(N,n)] by means of the shifted binomial transform

b_n = sum_{k=1}^{n} (-1)^(n-k) * nchoosek(n-1,k-1) * b_k^- for n >= 1 b_0 = 1

Input:

bm: vector of length n+1 holding b_0^-,…,b_n^-, i.e. bm(i) is the moment: of order i-1 and bm(1) = b_0^- = 1

Output:

b: vector of length n+1 holding b_0,…,b_n, with the same orientation: as bm

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, eq. (14).

Example: b = moment_binomial_from_negbinomial([1,2,4,22/3])

moment_binomial_from_factorial(f)

b = moment_binomial_from_factorial(f)

Converts the factorial moments f_n = E[N(N-1)…(N-n+1)] of a discrete random variable N into the binomial moments b_n = E[nchoosek(N,n)] via the one-to-one correspondence

b_n = f_n / n!

Input:

f: vector of length n+1 holding f_0,…,f_n, i.e. f(i) is the moment of: order i-1 and f(1) = f_0 = 1

Output:

b: vector of length n+1 holding b_0,…,b_n, with the same orientation: as f

Reference: A. Heindl and A. van de Liefvoort. Moment conversions for discrete distributions. PMCCS, 2003, Section 4.

Example: b = moment_binomial_from_factorial([1,2,4,8])