The Combinatorics of Distributions

Bernd Günther

Distributions, which are the various ways of distributing a certain number of objects of different classes among a collection of targets, have been the subject of combinatorial investigations since MacMahon’s 1917 monograph. In this paper we apply them to a simulation of superimposed random coding. Furthermore, asymptotic estimates are provided using logarithmic polynomials (related to the well-known Bell polynomials) for symbolic and numeric calculation.

Introduction

We denote by B(m,w,k) the number of ways m classes of w indistinguishable objects can be distributed among k different boxes, leaving no box empty and such that no box contains more than one object from the same class. In MacMahon’s terms [] (cf. also []) we speak of the distribution of objects of specification {(w^m)} into boxes of specification {(1^k)}.

Relaxing the requirement that no box should remain empty gives (((({{k}, {w}})))^m) distributions. Since each subcollection consisting of j boxes has B(m,w,j) permissible distributions, we obtain

(((({{k}, {w}})))^m)=(∑_(j=w)^k)(({{k}, {j}})) B(m,w,j).

Equation () can be readily solved:

B(m,w,k)=(∑_(j=w)^k)(((-1))^(k-j))(({{k}, {j}}))(((({{j}, {w}})))^m),

which holds for m>0 with B(1,w,k)=(Subscript[δ, w k]) for m=1.

numdistrib[m_Integer,w_Integer,k_Integer]:=(∑_(j=w)^k)(-1)^(k-j) Binomial[k,j] Binomial[j,w]^m /;m>0

It is convenient to extend the definition to the case m=0by setting B(0,w,k)=(Subscript[δ, k 0]), where, of course, we tacitly assume w>0.

numdistrib[0,w_Integer,k_Integer]:=KroneckerDelta[k,0]

For example, with m=2, w=3, k=4, there are 12 admissible distributions.

numdistrib[2,3,4]

Applications of distributions are abundant. Markov [] considered a lottery performing regular drawings of w out of n lots, and asked for the probability (Subscript[p, k]) that the number of lots that have shown up at least once after m drawings equals k. Here the lots take over the role of the boxes, which are occupied by an object of class j if it is drawn in the (j^th) round. Hence

(Subscript[p, k])=B(m,w,k)(({{n}, {k}}))(((({{n}, {w}})))^(-m)).

Superimposed Random Coding [, ] is an application to coding theory. Consider a collection of randomly generated binary code words of length n and weight w (i.e., the number of 1-bits equals w). A set of code words (Subscript[c, 1]),… ,(Subscript[c, m]) can be combined by superposition c=(Subscript[c, 1])∨ ⋯∨(Subscript[c, m]), so that c contains a 1-bit in all those places where at least one of the (Subscript[c, i]) has one. Notice that c has weight k with probability (Subscript[p, k]) as in (). For “decoding,” we tentatively assume that the (j^th) code word is present if its 1-bits are among those of c. This way we will certainly not miss any correct matches but may encounter a few false ones; the error probability must be investigated.

In practice we are testing for the presence of a set of code words (Subscript[c, 1])',…, (Subscript[c, s])' simultaneously by checking if its superposition c'=(Subscript[c, 1])'∨⋯∨(Subscript[c, s])' is covered by c. c' is of weight ℓ with probability

(Subscript[p, ℓ])'=B(s,w,ℓ)(({{n}, {ℓ}}))(((({{n}, {w}})))^(-s)).

Since there are (({{n}, {ℓ,k-ℓ,n-k}})) pairs of subsets of size ℓ and k (the first one contained in the second), the probability of a purely accidental match evaluates to

(Subscript[P, error])=(∑_(ℓ≤w))(({{n}, {ℓ,k-ℓ,n-k}}))(((({{n}, {ℓ}})))^(-1))(((({{n}, {k}})))^(-1))(Subscript[p, ℓ])'(Subscript[p, k]).

More details and a simulation are given in the supplemental notebook, which is available from www.mathematica-journal.com/data/uploads/2011/07/CombDistrib3Suppl.nb.

Asymptotics

Even for moderate values of m, w, and k, the number of distributions B(m,w,k) increases very quickly and soon exceeds the possibilities of ordinary 64-bit arithmetic. Fortunately, an asymptotic formula is available.

Theorem 1

If k,w→∞ such that the ratio (k/(m w))=t remains constant with 0<t<1 and if m either remains constant or also m→∞, but  m∈O((w^α)) for suitable α<(1/5), then

B(m, w, k) ∼((((m^(m w))(((1/(m(1-ζ)))-(1/m))((((1/(1-(ζ^m)))-1))^(-((1-(ζ^m))/(m(1-ζ)))))))^(m w))/(Sqrt[(((2 π w))^m)(ζ^m)(1-((1-(ζ^m))/(m(1-ζ)))(ζ^(m-1)))])),

where the new parameter 0<ζ<1 is defined by

((1-(ζ^m))/(m(1-ζ)))=(1/m) (∑_(j=0)^(m-1))(ζ^j)=t.

Sketch of Proof

The number of distributions can be computed by induction on m as

(Underoverscript[B(m+1,w,k)=∑, (j=k-w), k])(({{k}, {k-j,w-k+j,k-w}}))B(m,w,j),

which readily leads to

B(m,w,k)=(∑_(w=(Subscript[j, 1])≤(Subscript[j, 1])≤⋯≤(Subscript[j, m])=k))(∏_(i=2)^m)(({{(Subscript[j, i])}, {(Subscript[j, i])-(Subscript[j, i-1]),w-(Subscript[j, i])+(Subscript[j, i-1]),(Subscript[j, i])-w}})).

We resolve the multinomial coefficient as a product of factorials, apply Stirling’s formula, and approximate the sum by an integral over the variables (Subscript[ξ, i])=((Subscript[j, i])/w):

B(m,w,k)∼∫⋯(∫_((Subscript[ξ, 2])≤⋯≤(Subscript[ξ, m-1])))(Subscript[φ, w])((Subscript[ξ, 2]),⋯,(Subscript[ξ, m-1])) d(Subscript[ξ, 2])⋯d(Subscript[ξ, m-1]),

(Subscript[φ, w])((Subscript[ξ, 2]),⋯,(Subscript[ξ, m-1]))=(((2π w))^(-((m-1)/2)))(Sqrt[(∏_(i=2)^m)((Subscript[ξ, i])/(((Subscript[ξ, i])-(Subscript[ξ, i-1])) (1-(Subscript[ξ, i])+(Subscript[ξ, i-1])) ((Subscript[ξ, i])-1)))])(((ψ((Subscript[ξ, 2]),⋯,(Subscript[ξ, m-1]))))^w),

ψ((Subscript[ξ, 2]),⋯,(Subscript[ξ, m-1]))=(∏_(i=2)^m)(((Subscript[ξ, i])^(Subscript[ξ, i]))/(((((Subscript[ξ, i])-(Subscript[ξ, i-1])))^((Subscript[ξ, i])-(Subscript[ξ, i-1]))) (((1-(Subscript[ξ, i])+(Subscript[ξ, i-1])) )^(1-(Subscript[ξ, i])+(Subscript[ξ, i-1])))((((Subscript[ξ, i])-1))^((Subscript[ξ, i])-1)))).

The function ψ has a global maximum in the interior of the integration domain at (Subscript[ξ, i]^( * ))=((1-(ζ^i) )/(1-ζ)) with ζ chosen as in the theorem. Under such circumstances, (ψ^w) behaves like a Dirac delta distribution as w->∞:

(((ψ((Subscript[ξ, 2]),⋯,(Subscript[ξ, m-1]))))^w)∼(Subscript[β, w])δ((Subscript[ξ, 2])-(Subscript[ξ, 2]^( * )),…,(Subscript[ξ, m-1])-(Subscript[ξ, m-1]^( * ))),

(Subscript[β, w])=(Sqrt[(((((2 π))^(m-2))((ψ((Subscript[ξ, 2]^( * )),…,(Subscript[ξ, m-1]^( * ))))^(m+2 w-2)))/((w^(m-2))|(((∂^2)ψ)/(∂(Subscript[ξ, i])∂(Subscript[ξ, j])))|))]),

where |(((∂^2)ψ)/(∂(Subscript[ξ, i])∂(Subscript[ξ, j])))| denotes the absolute value of the Hessian of ψ at the location of the maximum. The validation of our estimates requires the interchange of limit operations at several instances, which can be justified if m∈O((w^α)). More details of the proof cannot be given within the present scope. We recommend [] as general reference for asymptotics.

Notice that the maximal possible value t=1 corresponding to k=m w, where every object is assigned to a different box, is not covered by the theorem.

We are interested in the limit case m→∞, where the number of classes used to typify our objects is very large; imagining that the different classes are distributed one after another (m can be interpreted as the time parameter of a Markov chain), this means that we want to anticipate the far future. However, the case m→∞, though included in the theorem, is not handled very conveniently, most notably because of the necessity of solving the (m-1)-degree polynomial equation (); some other deficiencies are described in the supplemental notebook.

Setting ϱ=m (1-ζ) and ε=(1/m) in equation () is equivalent to

log[1- t ϱ]=(1/ε)log[1-ε ϱ],

which is an analytic equation, since the singularity at ε=0 is inessential. By the implicit function theorem ϱ can be expressed as a power series in ε, and we are going to compute the coefficients. The logarithmic polynomials will be required as tools, for whose introduction we make a small digression.

Logarithmic Polynomials

We recall that the logarithmic polynomials (Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n])] are defined by

(Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n])]=(∑_(k=1)^n)(((-1))^(k-1)) (k-1)! (Subscript[B, n,k])[(Subscript[x, 1]),… (Subscript[x, n-k+1])],

where (Subscript[B, n,k]) denotes the Bell polynomials familiar from Faà di Bruno’s formula []. Observe that (Subscript[L, n]) is a polynomial in n indeterminates with integer coefficients. Their sequence satisfies the condition

log[1+(∑_(n=1)^∞)(((Subscript[x, n]) (t^n))/(n!))]=(∑_(n=1)^∞)(((t^n) (Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n])])/(n!)).

The following implementation is basically an adaption of [], exercise 16b (we are indebted to the referee).

logpoly[x_List]:= Module[{n,listOfPartitions,monomial,homdeg,aFactor},
n=Length[x];
listOfPartitions=Map[Tally,IntegerPartitions[n]];
aFactor[anEntry_]:=(x[[anEntry[[1]]]]/anEntry[[1]]!)^(anEntry[[2]])/anEntry[[2]]!;
monomial[aPartition_]:=(homdeg=Plus@@Map[Last,aPartition];
(-1)^(homdeg-1) n! (homdeg-1)! Times@@Map[aFactor,aPartition]);
Total@ Map[monomial,listOfPartitions]]/;Length[x]>0

Here are the first five log polynomials.

TableForm[Table[logpoly[Take[{(Subscript[x, 1]),(Subscript[x, 2]),(Subscript[x, 3]),(Subscript[x, 4]),(Subscript[x, 5])},n]],{n,5}],TableHeadings->Automatic]

Substituting (t^k) (Subscript[x, k]) for (Subscript[x, k]) in () one easily derives

(Subscript[L, n])[t (Subscript[x, 1]),(t^2)(Subscript[x, 2,])…,(t^n)(Subscript[x, n])]=(t^n)(Subscript[L, n])[(Subscript[x, 1]),…,(Subscript[x, n])],

thus the logarithmic polynomials are not homogeneous but of uniform weight, where the variable (Subscript[x, k]) counts with weight k. Furthermore, taking partial derivatives with respect to the highest argument,

(∂/(∂(Subscript[x, n]))) (Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n])]≡1,

which means that (Subscript[L, n]) depends additively on the last variable:

(Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n])]=(Subscript[L, n])[(Subscript[x, 1]),… (Subscript[x, n-1]),0] + (Subscript[x, n]).

Asymptotics, Resumed

Thus equipped, we set out to solve equation () for ϱ=(∑_(n=0)^∞)(((Subscript[ϱ, n])[t] (Subscript[(ε^n),  ]))/(n!)). The constant coefficient (Subscript[ϱ, 0])[t] can be obtained from () by taking limits:

log(1-t (Subscript[ϱ, 0]))=-(Subscript[ϱ, 0]).

It will be convenient to introduce η=exp(-(Subscript[ϱ, 0]))=1-t (Subscript[ϱ, 0]), 0<η<1, as a new independent parameter, so that

t=-((1-η)/(log η)).

This puts us in a position to put () into a form where () can be applied, namely where the constant coefficients of the involved power series equal 1:

log[1-(∑_(n=1)^∞)((t (Subscript[ϱ, n]) (Subscript[(ε^n),  ]))/(η n!))]=((1 )/ε)log[1-(∑_(n=1)^∞)(((Subscript[ϱ, n-1]) (Subscript[(ε^n),  ]))/((n-1)!))]+ (Subscript[ϱ, 0]).

By comparing coefficients we obtain

(Subscript[L, n])[-((t (Subscript[ϱ, 1]))/η),…,-((t (Subscript[ϱ, n]))/η)]=(1/(n+1)) (Subscript[L, n+1])[-(Subscript[ϱ, 0]),-2(Subscript[ϱ, 1]),…,-(n+1) (Subscript[ϱ, n])].

It will be seen presently that this determines (Subscript[ϱ, n]) as a rational expression in t and η; we observe that these two quantities, transcendentally related by (), are algebraically independent. To factor out the denominator, we set (Subscript[ϱ, n])=((η (((1-η))^(n+1)))/(((t^(n+1 ))(t-η))^(2n-1))) (Subscript[σ, n]) with (Subscript[σ, 0])=(1/(η (t-η))). Observing (), we obtain

(Subscript[L, n])[-(t-η) (1-η) (Subscript[σ, 1]),…,-(t-η) (1-η) (Subscript[σ, n])] = (((1-η))/((n+1) t (((t-η))^2))) (Subscript[L, n+1])[-(((t-η))^2) ,-2 η (((t-η))^(3 )) (Subscript[σ, 1]),…,-(n+1) η (((t-η))^(3 )) (Subscript[σ, n])].

Now from ():

(Subscript[σ, n])=(t/((1-η) (((t-η))^2)))(Subscript[L, n])[-(t-η) (1-η) (Subscript[σ, 1]),…,-(t-η) (1-η) (Subscript[σ, n-1]),0] 
-(1/((n+1) (((t-η))^4))) (Subscript[L, n+1])[-(((t-η))^2) ,-2 η (((t-η))^(3 )) (Subscript[σ, 1]),…,-n η (((t-η))^(3 )) (Subscript[σ, n-1]),0].

By uniformity of weight (), the denominators cancel, thus revealing (Subscript[σ, n]) for n≥1 as a polynomial.

sigmapoly[k_,η_,t_]:= sigmapoly[k,η,t]= Collect[Simplify[t logpoly[Append[Table[-(t-η) (1-η) sigmapoly[j,η,t],{j,1,k-1}],0]] / ((1-η) (t-η)^2) - logpoly[Append[Table[If[j==1,-(t-η)^2,-j η (t-η)^3 sigmapoly[j-1,η,t]],{j,1,k}],0]]/ ((k+1) (t-η)^4)],t]

Here are the first three cases of (Subscript[σ, n]).

TableForm[Table[sigmapoly[k,η,t],{k,3}],TableHeadings->Automatic]

We can now check the correctness of our coefficients by inserting them into equation ().

maxDim=5;
rhotest=(1-η)/t +Sum[(η*(1 - η)^(k + 1) * sigmapoly[k,η,t] *ε^k)/(k! t^(k + 1)*(t - η)^(2*k - 1)) ,{k,maxDim}];Simplify[Log[1-t rhotest +O[ε]^(maxDim+1)] - Log[1-ε rhotest + O[ε]^(maxDim+2)]/ε]

This is 0 because of (). Let us summarize the results achieved so far.

Lemma 1

The unique solution ζ>0 of the equation ((1-(ζ^m))/(m(1-ζ)))=t is given by the power series

ζ=1-((1-η)/(m t))-(∑_(n=1)^∞)((η (((1-η))^(n+1)) (Subscript[σ, n]))/(n! (((t^(n+1 ))(t-η))^(2n-1)) (m^(n+1)))),

The rest now is fairly standard. We observe that in (),

log(((1/(m(1-ζ)))-(1/m))((((1/(1-(ζ^m)))-1))^(-((1-(ζ^m))/(m(1-ζ))))))=log((1/ϱ)-ε)-t log((1/(t ϱ))-1]),

and we can easily compute a power series expansion of the expression log((1/ϱ)-ε)-t log((1/(t ϱ))-1)=(∑_(n=0)^∞)(Subscript[α, n])(ε^n) from the one of ϱ.

maxDim=5;
rhotest=(1-η)/t +Sum[(η*(1 - η)^(k + 1) * sigmapoly[k,η,t] *ε^k)/(k! t^(k + 1)*(t - η)^(2*k - 1)) ,{k,maxDim}];pbase=Simplify[Log[1/rhotest -ε+O[ε]^(maxDim+1)] - t Log[1/(t rhotest)-1 + O[ε]^(maxDim+1)]];
α[n_]:=Coefficient[pbase,ε,n];

The first two coefficients are as follows.

Exp[α[0]]/.t->-(1-η)/Log[η]

Exp[α[1]]/.t->-(1-η)/Log[η]

This provides us with complete control over the numerator of (). To handle the denominator we observe (ζ^m)=1-t m (1-ζ) and therefore, from (), (lim_(m->∞))(ζ^m)=η. We also observe (lim_(m->∞))ζ=1, because ζ depends implicitly on m. In consequence,

(lim_(m->∞))(Sqrt[(ζ^m)(1-((1-(ζ^m))/(m(1-ζ)))(ζ^(m-1)))])=(Sqrt[η(1-t η)]).

In using this limit approximation we accept a slight penalty of accuracy in favor of a simpler formula. Assembling all parts, we have proved the following theorem.

Theorem 2

If k,m,w→∞ such that the ratio (k/(m w))=t remains constant with 0<t<1 and m∈O((w^α)) for suitable α<(1/5), then

B(m, w, k)∼(Sqrt[((η^(w-1))/((((2 π w))^m) (1-t η)))]) (((m ((((((1-η)/η)))^((1-η)/(-log η)))/(-log η))))^(m w))exp((∑_(n=2)^∞)(((Subscript[α, n])w)/(m^(n-1)))),

Since m is not supposed to grow faster than (w^(1/5)), at least the first five terms of the power series in () will be significant, but higher terms may be negligible depending on the actual growth rate.

The derivation of error bounds for the asymptotic relation () is beyond our present scope. A numerical example showing that the approximation is rather good is contained in the supplemental notebook.

References

P. A. MacMahon, Combinatory Analysis, 3rd ed. reprint, Providence, RI: AMS Chelsea, 2001.

J. Riordan, Introduction to Combinatorial Analysis, Mineola, NY: Dover Publications, 2002.

A. A. Markov, Wahrscheinlichkeitsrechnung, Leipzig, Berlin: Teubner, 1912. catalog.hathitrust.org/Record/012107751.

B. Günther, “On the Probability Distribution of Superimposed Random Codes,” IEEE Transactions on Information Theory, 54(7), 2008 pp. 3206–3210. doi:10.1109/TIT.2008.924658.

C. S. Roberts, “Partial-Match Retrieval via the Method of Superimposed Codes,” Proceedings of the IEEE, 67(12), 1979 pp. 1624–1642. ieeexplore.ieee.org/xpl/freeabs_all.jsp?arnumber=1455812.

N. G. de Bruijn, Asymptotic Methods in Analysis, Amsterdam: North-Holland Publishing Co.; New York: Interscience Publishers, 1958.

L. Comtet, Advanced Combinatorics—The Art of Finite and Infinite Expansions, rev. and enl. ed. (J. W. Nienhuys, trans.), Dordrecht, Boston: D. Reidel Publishing Company, 1974 p. 140.

M. Trott, The Mathematica GuideBook for Numerics, Berlin: Springer, 2005.