Grouping lists together in a proportional election: image of a Dirichlet distribution by the D'Hondt method of proportional allotment

Question

A real-world motivation for this question is given below. But first let me recall what the Jefferson-D'Hondt “greatest divisors” method of proportional allotment (often used in electoral systems to elect an assembly) means in mathematical terms:

The D'Hondt method: Given a (sufficiently general) point $P$ on the $(n-1)$-dimensional simplex, viꝫ. $P = (p_1,\ldots,p_n)$ with $p_i$ nonnegative and summing to $1$ (the “votes”), and given $N \geq 1$ (the number of “seats” to allot), we define a point $\tau_N(P)$ on the “$N$-discrete $(n-1)$-dimensional simplex”, that is, an $n$-tuple $(m_1,\ldots,m_n)$ of natural numbers with sum $N$, as follows: $\tau_N(P)$ is the $n$-tuple $(m_1,\ldots,m_n)$ such that the $m_i = \lfloor t p_i\rfloor$ have sum exactly $N$, for some positive real $t$. In more geometric terms¹, consider the half-ray $\{tP : t>0\}$ through $P$ in $\mathbb{R}^n$ and see which cube $\prod_{i=1}^n [m_i,m_i+1]$ with $m_1+\cdots+m_n=N$ it intersects (even more visually: we are looking at the arrangement of unit cubes corresponding to the $N$-discrete simplex, and the division of the simplex given by $\tau_N$ corresponds to looking at this arrangement from the origin). See the Wikipedia article linked above for more computational descriptions. This $\tau_N(P)$ is well-defined except for cases of “ties” in the votes, which are contained in a finite set of hyperplanes in the simplex, and since I will be considering probabilities in my question, these are irrelevant here.

The D'Hondt method has the following interesting property, for which I do not know a reference in the literature (but here is a proof that was given to me long ago on sci.math.research): if $P$ is drawn uniformly on the simplex, then $\tau_N(P)$ follows a uniform distribution on the $N$-dimensional simplex. (In slightly handwavy terms: “if you don't know anything about the result of the vote, you don't know anything about the attribution of the seats”.) This is what leads me to hope that the following question might not be too untrackable:

General question: Suppose $P$ is drawn from a Dirichlet distribution on the simplex, what distribution does $\tau_N(P)$ follow on the $N$-discrete simplex? What are the expected values of its coordinates?

Now since this question might turn out to be too difficult for explicit computations, here is a simpler form that I hope is made more manageable by considering Dirichlet distributions where all parameters are $1$ except one which is integer:

Hopefully simpler version: Let $1\leq k\leq n$, and consider the map $\pi\colon (p_1,\ldots,p_n) \mapsto (p_1,\ldots,p_{n-k}, p_{n-k+1}+\cdots+p_n)$ from the $(n-1)$-dimensional simplex to the $(n-k)$-dimensional simplex (which sums the last $k$ coordinates), so that $\pi(P)$ follows a Dirichlet distribution with parameters $(1,\ldots,1,k)$ if $P$ is uniform: what is the expected value of (the last coordinate of) $\tau_N(\pi(P))$ if $P$ is uniform?

Note that if we call $\pi_N$ the analogous map $(m_1,\ldots,m_n) \mapsto (m_1,\ldots,m_{n-k}, m_{n-k+1}+\cdots+m_n)$ on the $N$-discrete simplices, then $\pi_N(\tau(P))$ has expected value $(\frac{N}{n},\ldots,\frac{N}{n},\frac{kN}{n})$ since $\tau(P)$ is uniform by the property remarked above. The expected value of the last coordinate of $\tau_N(\pi(P))$ is greater than $\frac{kN}{n}$ and the question is “by how much?”.

Real-life motivation: European parliament elections are happening in 2024, and France allots its seats by the D'Hondt method. One of the questions that come up in the political debate is how many seats a certain set of $k$ political parties might gain if they run under a single list instead of $k$ separate lists.

Experimental results: By sampling points on the unit simplex and guessing the denominators, I have arrived at the following results for expected value of the last coordinate of $\tau_N(\pi(P))$ minus $\frac{kN}{n}$:

• for $(n,k)=(3,2)$ we find $\frac{N}{6(N+1)}$,

• for $(n,k)=(4,3)$ we find $\frac{N}{4(N+1)}$,

• for $(n,k)=(4,2)$ we find $\frac{N(3N+5)}{12(N+1)(N+2)}$,

• for $(n,k)=(5,4)$ we find $\frac{N(9N^2+17N+7)}{30(N+1)^3}$,

• for $(n,k)=(5,3)$ we find $\frac{N(24N^3+114N^2+177N+89)}{60(N+1)^2(N+2)^2}$,

• for $(n,k)=(5,2)$ we find $\frac{N(3N^2+13N+13)}{10(N+1)(N+2)(N+3)}$

(in other words, the quantity above is the expected advantage that $k$ party lists out of $n$ have in grouping together if there are $N$ seats to fill and the result of the election is random).

Integral expression: [added 2023-05-11] I should mention that the distribution of $\tau_N(P)$ in answer to the “general question” is proportional to the following integral expression (for $(m_1,\ldots,m_n)$ an $n$-tuple of natural numbers with sum $N$): $$ \sum_{j=1}^n m_j \int_{\prod_{i\neq j}[m_i,m_i+1]} \frac{x_1^{\alpha_1-1}\cdots x_n^{\alpha_n-1}}{(x_1+\cdots+x_n)^{\alpha_1+\cdots+\alpha_n}} dx_1\cdots dx_n $$ (a straightforward but tedious translation of the description of the D'Hondt method as projecting cubes to the simplex: the sum over $j$ ranges over the facets of the cube $\prod_{i=1}^n[m_i,m_{i+1}]$); here, $x_j$ in the integral should be interpreted as being identically $m_j$, and $dx_j$ should be omitted. This is the integral which has been used (laboriously) to compute the “experimental results” indicated above.

1. Yet another description consists of noting that (again, for sufficiently general $P$), the ray $\{tP\}$ cuts the walls of the subdivision of space into unit cubes in a pattern that forms a Sturmian sequence, and the D'Hondt method counts the letters of each kind within the $N$ first letters in this sequence.

Answer 1

Here is a way. Start from the exact form of your integral expression (i.e. supply the "Dirichlet factor" $\frac{\Gamma(\alpha_1+\ldots+\alpha_n)}{\Gamma(\alpha_1)\cdots\Gamma(\alpha_n)})$.

(1) Writing $\frac{1}{(x_1+\ldots+x_n)^{\alpha_1+\ldots+\alpha_n}}=\int_0^{\infty}\frac{ t^{\alpha_1+\ldots+\alpha_n-1}}{\Gamma(\alpha_1+\ldots+\alpha_n)}\,e^{-t(x_1+\ldots+x_n)}\,dt$ and changing the order of integration (and denoting by $p_{n,N}$ the distribution of $\tau(N)(P)=:(M_{1,N},\ldots,M_{n,N})$) leads to $$p_{n,N}(m_1,\ldots,m_n)=\sum_{j=1}^n \int_0^\infty \frac{m_j^{\alpha_j}t^{\alpha_j-1}}{\Gamma(\alpha_j)}\prod_{i\neq j}\big (R_{\alpha_i}(m_it)-R_{\alpha_i}((m_i+1)t)\big) e^{-m_jt}\,dt$$ where $R_{\alpha} (t):=\int_t^\infty \frac{y^{\alpha-1}}{\Gamma(\alpha)} e^{-y}\,dy$ (and $m_1+\ldots+m_n=N$).

(2) For $\alpha=k$ (a positive integer) $$R_k(t)=q_k(t) e^{-t} \mbox{ with } q_k(t)=\sum_{i=0}^{k-1}\frac{t^i}{i!}\;\;.$$ If all $\alpha_i$ are all integers the integrals above therefore lead to explicit rational expressions for the probabilities $p(m_1,\ldots,m_n)$. For the case $\alpha_1=\ldots=\alpha_{n-1}=1, \alpha_n=k$ (positive integer) one gets \begin{align*} p_{n,N}(m_1,\ldots,m_n)= &\frac{m_n^k}{(k-1)!}\int_0^\infty t^{k-1} e^{-Nt} (1-e^{-t})^{n-1}\,dt\\& + (N-m_n)\int_0^\infty \big(q_k(m_nt)-e^{-t}q_k((m_n+1)t)\big) e^{-Nt}(1-e^{-t})^{n-2}\,dt\end{align*} ADDED: a fully explicit solution for the case $\alpha_1=\ldots=\alpha_{n-1}=1, \alpha_n=k$ (positive integer).

(3) In the sequel $n\geq 2$, and $k$ is fixed and suppressed from the notation. Working out the integrals above gives (for any nonnegative integers $m_1,\ldots,m_{n-1},m$ with $m_1+\ldots +m_n+m=N$) \begin{align*} p_{n,N}(m_1,\ldots,m_{n-1},m)=I_{n-1,N}(m) + J_{n-2,N}(m) \end{align*} where \begin{align*} I_{n,N}(m)&=m^k\sum_{j=0}^{n}(-1)^j {n \choose j}\frac{1}{(N+j)^k}\\ J_{n,N}(m)&=\sum_{j=0}^{n} (-1)^j {n \choose j} \frac{N-m}{N-m+j}\Big(\big(\frac{m+1}{N+j+1}\big)^k-\big(\frac{m}{N+j}\big)^k\Big)\\ \end{align*}

Since the joint distribution $p_{n,N}(m_1,\ldots,m_{n-1},m)$ depends on $m_1,\ldots,m_{n-1}$ only through $m_1+\ldots+m_{n-1}=N-m$ (we will from now on simply write $p_{n,N}(m)$, and) the distribution $f_{n,N}(m):=\mathbb{P}(M_{n,N}=m)$ of the last coordinate $M_{n,N}$ is simply given by \begin{align*} f_{n,N}(m)= {N-m+n-2 \choose n-2} \cdot p_{n,N}(m) \end{align*} (since $ {N-m+n-2 \choose n-2}$ is the number of $n-1$-tuples of nonnegative integers summing to $N-M$).

(4) Examples

(a) For $k=1$ both integrals above are the same, and we get \begin{align*} p_{n,N}(m)=N\int_0^\infty e^{-Nt} (1-e^{-t})^{n-1}\,dt=N\, B(N-1,n-1)=\frac{N!\,(n-1)!}{(N+n-1)!}\;\;, \end{align*} confirming that $\tau(N)(P)$ is uniformly distributed on the discrete simplex (as was already noted above).

(b) For $n=2$ we find $I_{1,N}=m^k\big(\frac{1}{N^k}-\frac{1}{(N+1)^k}\big), J_{0,N}=(\frac{m+1}{N+1})^k-(\frac{m}{N})^k$ so that

\begin{align*} f_{2,N}(m)=\frac{(m+1)^k-m^k}{(N+1)^k}=p_{2,N}(m),\;\;\;\;m=0,\ldots,N \end{align*}

and similar computations give for $n=3$

\begin{align*} f_{3,N}(m)=& \Big((m+1)^k-m^k\Big)\Big(\frac{N+2}{(N+1)^k}-\frac{N+1}{(N+2)^k}\Big) \\&-\Big((m+1)^{k+1}+m^{k+1}\Big)\Big(\frac{1}{(N+1)^k}-\frac{1}{(N+2)^k}\Big) \end{align*} (5) The distribution of $M_{n,N}$ for arbritrary $n$. Although the above description of the distribution of $\tau(N)(P)$ resp. $M_{n,N}$ is (in principle) complete, it is not very transparent. A surprisingly simple description can be obtained if one uses finite differences. Let \begin{align*} \delta_n(N):=\sum_{j=0}^n{n \choose j} \frac{(-1)^j}{(N+j)^k}\end{align*} (so that $(-1)^n \delta_n(N)$ is the $n$-th forward difference of $1/x^k$ at $N$). We then have

Lemma \begin{align*} (1)&\;\;{N-m+n \choose n}J(n,N)=\sum_{j=0}^n {N-m-1+j \choose j}\Big((m+1)^k\,\delta_j(N+1) -m^k \delta_j(N)\Big)\;\;\\ (2)&\;\;f_{n+2,N}(m)=\sum_{j=0}^n\Big[(m+1)^k{N-(m+1)+j \choose j}-m^k{N-m+j \choose j}\Big]\;\delta_j(N+1)\;\;\;. \end{align*} Proof (Sketch) (1): use (show) the equality $$\sum_{j=i}^n {N-m-1+j \choose j}{ j \choose i}=\frac{N-m}{N-m+i}{N-m+n \choose n} {n \choose i}\;\;\;.$$ (2): the coefficients of $(m+1)^k$ follow directly from (1). The coefficients of $m^k$ follow by using $I(n,N)=m^k\delta_{n+1}(N)$, $\delta_{n}(N)=\delta_{n-1}(N)-\delta_{n-1}(N+1)$ and the addition rule for binomial cofficients.;- End Proof

(6) The expectation of $M_{n,N}$

From (2) of the lemma $\mathbb{P}(M_{n+2,N}\leq m)=\sum_{j=0}^n (m+1)^k{N-(m+1)+j \choose j}\;\delta_j(N+1)$ and thus \begin{align*} \mathbb{E}(M_{n+2,N})&=\sum_{m=0}^{N-1} \mathbb{P}(M_{n+2,N}>m)\\ &=\sum_{m=0}^{N-1} \big(1-\mathbb{P}(M_{n+2,N}\leq m)\big)\\ &=N -\sum_{j=0}^n\Big[\sum_{m=1}^N m^k{N-m+j \choose j}\Big]\;\delta_j(N+1) \end{align*} The sums $\sum_{m=1}^N m^k{N-m+j \choose j}$ can be expressed in various ways.

Denoting the Stirling numbers of the second kind by ${k \brace d}$, the falling factorials by $(m)_d$ and using $m^k=\sum_{d=1}^k { k \brace d} (m)_d$ and $\sum_{m=1}^N (m)_d {N-m+j \choose j}=d!{N+j+1 \choose d+j+1}$ gives \begin{align*} \sum_{m=1}^N m^k {N-m+j \choose j}=\sum_{d=1}^k d!\,{ k \brace d}{N+j+1 \choose d+j+1} \end{align*} Plugging this in above gives an explicit formula for the expectation. Alternatively one could could expand the binomial cofficient as a polynomial in $m$, and use power sums/Faulhabers formula (as has been done in the solution below).

Answer 2

With another approach we get

$(n, k) = (k+1, k)\;\implies\;\mathrm{E}[m_0] - \frac{kN}{n} = $ $$\frac{N}{k+1} - \frac{1}{(1+N)^k}\sum_{i=1}^N i^k$$ $(n, k) = (k+2, k)\;\implies\;\mathrm{E}[m_0] - \frac{kN}{n} = $ $$\frac{2N}{k+2} - \left(\frac{N+2}{(N+1)^k}-\frac{N+1}{(N+2)^k}\right)\sum_{i=1}^N i^k + \left(\frac{1}{(N+1)^k}-\frac{1}{(N+2)^k}\right) \sum_{i=1}^N i^{k+1} $$ where $m_0$ denotes your $m_{n-k+1}+\dots+m_n\,$. Applying Faulhaber's formula yields your five first results.

---

This is because $-$ ties put aside $-$ the numbers of seats $(m_0,\ldots,m_d)\in\mathbb{N}^{d+1}$ alloted to $d+1$ parties by the Jefferson-D'Hondt method is known to occur when the proportions of votes $(x_0, \dots, x_d)\in\mathbb{R}^{d+1}$ satisfy

$$\frac{m_i}{m_j+1}x_j \le x_i\qquad\forall i\ne j=0,\dots,d\;.$$

Choosing $x_d$ to stand for $1-x_0-\dots-x_{d-1}\,,\;$ these $d(d+1)$ inequalities upon $x_0,$ $\dots,$ $x_{d-1}$ intersect half-spaces of $\mathbb{R}^d$ and draw a rational bounded convex polytope $\mathcal{P} = \mathcal{P}(m_0,\dots,m_d)$ of $\mathbb{R}^d\,$.

The probability for $(m_0,\dots,m_d)$ to occur when the $x_i$'s follow a Dirichlet distribution with parameters $(\alpha_0,\dots,\alpha_d)$ is then $\mathrm{Pr}[m_0,\dots,m_d]=$

$$\frac{\Gamma(\alpha_0+\dots+\alpha_d)}{\Gamma(\alpha_0)\dots\Gamma(\alpha_d)}\int_{\mathcal{P}} x_0^{\alpha_0-1}\cdots x_{d-1}^{\alpha_{d-1}-1}(1-x_0-\dots-x_{d-1})^{\alpha_d-1} \,\mathrm{d}x_0\cdots \mathrm{d}x_{d-1}\;. $$

---

When $d=1$, writing $(x,y)$ or $(x,1-x)$ for $(x_0,x_1)\,$, and $(a,b)$ or $(a,N-a)$ for $(m_0,m_1)\,$, this is $$\frac{a}{b+1}(1-x)\le x\quad\text{and}\quad\frac{b}{a+1}x\le 1-x\;.$$ In other words $$x\in\mathcal{P} = \mathcal{P}(a,b) = \left[\frac{a}{N+1}, \frac{a+1}{N+1}\right]$$ so that $$ \mathrm{Pr}[a,b]= \frac{\Gamma(\alpha_0+\alpha_1)}{\Gamma(\alpha_0)\Gamma(\alpha_1)}\int_\frac{a}{N+1}^\frac{a+1}{N+1} x^{\alpha_0-1}(1-x)^{\alpha_1-1}\,\mathrm{d}x\;. $$

In the simpler version $n=1+k$ and $(\alpha_0,\alpha_1)=(k,1)$, we get $$ \mathrm{Pr}[a,b]= k\int_\frac{a}{N+1}^\frac{a+1}{N+1} x^{k-1}\,\mathrm{d}x = \left(\frac{a+1}{N+1}\right)^k - \left(\frac{a}{N+1}\right)^k $$ so that $$\begin{align*} \mathrm{E}[m_0]& = \sum_{a+b=N}a\,\mathrm{Pr}[a,b] = \frac{1}{(N+1)^k}\sum_{a=0}^N a\,\left((a+1)^k-a^k\right) = N - \frac{1}{(N+1)^k}\sum_{i=1}^N i^k \end{align*}$$

---

When $d=2$, indifferently writing $(x,y,z)$ or $(x,y,1-x-y)$ for $(x_0,x_1,x_2)\,$, and $(a,b,c)$ or $(a,b,N-a-b)$ for $(m_0,m_1,m_2)\,$, this is $$ \begin{array}{l} &\dfrac{a}{1+b}y\le x&\dfrac{a}{1+c}(1-x-y)\le x\\ \dfrac{b}{1+a}x\le y&&\dfrac{b}{1+c}(1-x-y)\le y\\ \dfrac{c}{1+a}x\le 1-x-y&\dfrac{c}{1+b}y\le 1-x-y& \end{array} $$

Here is a picture of the resulting $2$-dimensional $\mathcal{P}(a,b,c)$ polytopes when $N=a+b+c=6$: The vertices of $\mathcal{P}(a,b,c)$ are related to the following six binary strings and points of $\mathbb{R}^3$ $$\mathsf{100}\to\left(\tfrac{a+1}{N+1},\tfrac{b}{N+1},\tfrac{c}{N+1}\right),\\ \mathsf{101}\to\left(\tfrac{a+1}{N+2},\tfrac{b}{N+2},\tfrac{c+1}{N+2}\right),\qquad\mathsf{110}\to\left(\tfrac{a+1}{N+2},\tfrac{b+1}{N+2},\tfrac{c}{N+2}\right),\\ \mathsf{001}\to\left(\tfrac{a}{N+1},\tfrac{b}{N+1},\tfrac{c+1}{N+1}\right),\qquad\mathsf{010}\to\left(\tfrac{a}{N+1},\tfrac{b+1}{N+1},\tfrac{c}{N+1}\right),\\ \mathsf{011}\to\left(\tfrac{a}{N+2},\tfrac{b+1}{N+2},\tfrac{c+1}{N+2}\right).$$ They may not all be extreme points but they are always on the boundary. Edges are between any two points that happen to be at Hamming distance $1$.

In the simpler version $n=2+k$ and $(\alpha_0,\alpha_1,\alpha_2)=(k,1,1)$, the density is constant at $x$ constant and the length of the slice through the polytope $\mathcal{P}$ for $x$ constant is a piecewise linear function of $x$. Working out those linear terms yields the elementary integration $$\begin{align} \mathrm{Pr}[a,b,c] = k(k+1)\int_{\frac{a}{N+2}}^{\frac{a}{N+1}}x^{k-1}\left(\tfrac{N+2}{a}x-1\right)\mathrm{d}x + k(k+1)\int_{\frac{a}{N+1}}^{\frac{a+1}{N+2}}x^{k-1}\tfrac{1}{1+N-a}(1-x)\mathrm{d}x\qquad\\ + k(k+1)\int_{\frac{a+1}{N+2}}^{\frac{a+1}{N+1}}x^{k-1}\left(1-\tfrac{N+1}{a+1}x\right)\mathrm{d}x\\ = \left(\frac{a+1}{N+1}\right)^k - \frac{N-a}{N-a+1}\left(\frac{a+1}{N+2}\right)^k - \frac{N-a+2}{N-a+1}\left(\frac{a}{N+1}\right)^k + \left(\frac{a}{N+2}\right)^k \end{align}$$ which doesn't depend on $b$ or $c$ individually, just like the density of $(x,y,z)$ doesn't depend on $y$ or $z$ individually. For $N$ fixed we then get $$E[m_0] = \sum_{a+b+c=N} a\,\mathrm{Pr}[a,b,c] = \sum_{a=0}^Na\,(N-a+1)\,\mathrm{Pr}[a,0,N-a]$$ which yields the formula given in introduction.

---

I didn't figure out $d=3$ yet while experimenting with polymake. It is hard to come up with a description without many subcases, such as whether each $m_i$ is smaller or greater than $\frac{N+1}{2}$.