# Entropy of fair but correlated coin flips

Consider the joint distribution, $$p(\xi_1,...\xi_N)$$, with components defined as $$\xi_i=\mathrm{sign}(x_i)$$, with $$(x_1,...,x_N)\sim\mathcal{N}(0,\Sigma)$$ with $$\Sigma_{ij}=\delta_{ij}+(1-\delta_{ij})\tilde{\rho}$$, i.e. all off-diagonal entries of $$\Sigma$$ are $$\tilde{\rho}$$. Since all components $$x_i$$ have unit variance (i.e. diagonal entries of $$\Sigma$$ are $$1$$), $$\tilde{\rho}$$ is then the correlation of any pair $$(x_i,x_j)$$. Note that

1. all single component marginals $$p(\xi_i)=1/2$$, i.e. the coins are unbiased;
2. we can always set the value of $$\tilde{\rho}$$ so that any pair $$(\xi_i,\xi_j)$$ has the desired correlation of $$\rho$$.

What is the formula for the entropy of $$p(\xi_1,...\xi_N)$$, denoted $$H_N(\rho)$$?

The parametrization of $$\tilde{\rho}$$ by $$\rho$$ depends on $$N$$ and is obtained by calculating the pair marginal probability $$p(\xi_i=1,\xi_j=1)$$, denoted $$p_{11}$$. Let $$f_N(\tilde{\rho})$$ be the function of $$\tilde{\rho}$$ that gives this probability. The formula for the correlation between two binary random variables, $$\rho=\frac{p_{11}-p_{-1}p_1}{\sqrt{p_{-1}(1-p_{-1})p_1(1-p_1)}}={4}p_{11}-1\;,$$ with $$p_{\pm1}$$ denoting $$p(\xi_i=\pm 1)=p(\xi_j=\pm 1)=1/2$$, gives the desired parametrization, $$\tilde{\rho}_N(\rho):=f^{-1}_N\left(\frac{\rho+1}{4}\right)\;.\;\;\;\;(\mathrm{eq}.1)$$

Solution for $$N=2$$:

Computing the 2D Gaussian integral using spherical coordinates gives $$\tilde{\rho}_{N=2}(\rho)=\sin(\frac{\pi}{2}\rho)$$. In this special case, $$(\mathrm{eq}.1)$$ and symmetry constraints completely specify the distribution $$p(\xi_1,\xi_2)=(1+\xi_1\xi_2\rho)/4$$. Calculation of the entropy from its definition gives

$$H_{N=2}(\rho)=H(\xi_1,\xi_2)=\sum_{\eta=(1\pm\rho)/4}2\left(-\eta\log_2\eta\right).$$

We have $$H_{N=2}(0)=2$$ bits and $$H_{N=2}(1)=1$$ bit. In fact, we know $$H_N(\rho=0)=N$$ and $$H_N(\rho=1)=1$$ for all $$N$$.

For $$N>2$$:

There are lots of Gaussian integrals to do. Due to the symmetry in the index permutations there are only $$N$$ distinct integral values among the $$2^N$$ terms in the entropy. They can be grouped by how many 1s they contain. The pair of all $$-1$$s can be grouped with the singleton group of all $$1$$s due to the reflective symmetry in the plane normal to the main diagonal. We just need compute the multiplicity and the integral for each of the $$N$$ terms.

E.g. for $$N=2$$ there are 2 terms (listed above with multiplicity 2). For $$N=3$$ there are three terms (two for tuples containing one and two 1s), the third being the extreme group having $$(-1,-1,-1)$$ and $$(1,1,1)$$ in binary notation.

In section 3.2 of this paper, Six gives a recursive solution to the distribution. This could be used to compute the entropy and the parametrization function $$\tilde{\rho}_N(\rho)$$.

Accepted Answer: In the end, an alternative approach based on expressing $$x_i=\sqrt{1-\tilde{\rho}}y_i+\sqrt{\tilde{\rho}}s$$, with $$y_i,s\sim\mathcal{N}(0,1)$$, seems to be more straightforward and is the accepted answer below. That answer also indicates that the $$f_N(\tilde{\rho})$$ does NOT seem to depend on $$N$$ after all.

• Maybe I don't understand something, but is this really enough information for an answer? If $N=3$ and $\rho=0$ for example, we could have $\xi$ drawn from $\{0,1\}^3$ uniformly at random and the entropy is $3$, or we could have it drawn from the subset of tuples with an even number of $1$s uniformly at random and then the entropy is $2$.
– anon
Commented Apr 12, 2022 at 0:41
• The first example is how I am thinking about it. What is the motivation/reasoning for the 2nd example? I am confused by it. How are "even number of 1" tuples relevant to this problem? What about the odd ones? Order does matters here. Will try to clarify language. Commented Apr 12, 2022 at 1:19
• I showed how the information given (i.e. $N$ and $\rho$) is not enough to solve for the entropy, by giving two examples where $N$ and $\rho$ are the same but the entropies are different. The second example I gave is standard for showing how pairwise independence does not imply joint independence. In the second example, tuples with an odd number of $1$s have probability zero.
– anon
Commented Apr 12, 2022 at 1:32
• Thanks. Do you have a reference for such an example. I still haven't understood it (e.g. isn't (1,0,0) in the set of tuples and have an odd number of 1s?). I also added a note on how to generate samples from this distribution by sampling from multivariate normal. Commented Apr 12, 2022 at 1:48
• $(1,0,0)$ is in the set of all eight tuples, but it is not in the subset of four tuples which have an even number of $1$s. Why? Paul, Andre, and izimath provide this example in [this thread](math.stackexchange.com/questions/1783225/…), for instance. What haven't you understood about it? The sample space has literally only four outcomes each with probability one-fourth - it is straightforward to check that in this situation, each $\xi_i$ is a fair coin toss and each pair of them are independent (so $\rho=0$).
– anon
Commented Apr 12, 2022 at 1:57

Of course, the assumption $$H(\xi_i|\xi_{i+1},...,\xi_N)=H(\xi_i|\xi_{i+1})$$ seems wrong, even as an approximation. I doubt that there is a simple solution.

Here's at least some numerical results, from simulations. The horizontal axis is $$\rho$$, the vertical axis is the estimated entropy, normalized to $$[0,1]$$: $$H'=(H-1)/(N-1)$$

As you mention, the evaluation of the entropy amounts to integrate a multimensional gaussian (albeit zero mean and symmetric wrt the components) over all the quadrants. The answers to this question have some pointers. Eg this one

The most specific reference I found is this one .

Another interesting one: https://arxiv.org/abs/1009.2855

Added: from the same simulations, here are the conditional probabilities $$H(\xi_i|\xi_{i-1},...,\xi_1)$$ for diferent values of $$\rho$$. I wonder if there is non-zero limit for $$N\to \infty$$...

In case someone is interested, here's the Python code

import numpy as np
import math

def hb(p):
return -p * math.log2(p) - (1 - p) * math.log2(1 - p) if p > 0 and p < 1 else 0

def calca(nn, r):
return (math.sqrt((1 - nn) * r**2 + (nn - 2) * r + 1) - r + 1) / r

def calcrho(nn, aa):
return (nn + 2 * aa) / ((aa + 1) ** 2 + nn - 1)

def gen1(nn, aa):
y = np.random.normal(0, 1, nn)
return (np.sum(y) * np.ones(nn) + aa * y) / math.sqrt((aa + 1) ** 2 + (nn - 1))

def tryn(nn, rho, t):
aa = calca(nn, rho)
Pk = np.zeros((nn + 1, nn + 1))
PJk = np.zeros((nn + 1, nn + 1))  # joint prob
PCk = np.zeros((nn + 1, nn + 1))  # cond prob
for _ in range(t):
y = gen1(nn, aa)
z = np.array([1 if t > 0 else 0 for t in y])
c = 0
for k in range(nn):
c += z[k - 1] if k > 0 else 0
Pk[k][c] += 1
if z[k] == 1:
PJk[k][c] += 1
# normalize probabilities
for k in range(nn):
Pk[k] /= t
PJk[k] /= t
for k in range(nn):
for c in range(k + 1):
PCk[k][c] = PJk[k][c] / Pk[k][c] if Pk[k][c] > 0 else 0
r = calcrho(nn, aa)
print(f"{nn=} {r=} {aa=} {t=}")
# print(np.round(Pk[:5,:5],5))
# print(np.round(PJk[:5,:5],5))
# print(np.round(PCk,5))
H = np.zeros(nn)
for k in range(nn):
H[k] = 0
for c in range(k + 1):
H[k] += hb(PCk[k][c]) * Pk[k][c]
print('h', np.round(H, 5))

rho = 0.7
nn = 12
tryn(nn, rho, 800000)

• cool, thanks. So each goes to $1/N$ at $\rho=1$. And you sample from the Gaussian with the $sin(\pi/2\rho)$ covariance and then compute the entropy directly using empirical frequencies? Commented Apr 13, 2022 at 4:13
• Looks quadratic around $\rho=0$ for all $N$, but becomes more linear away from $\rho=0$ as $N$ is increased. It might be easier to see the large-$N$ limit if you subtract 1 and divide by $N-1$ to pin the $\rho=0,1$ values for all $N$ at $1$ and $0$ respectively. Commented Apr 13, 2022 at 4:23
• @puelmato Good point, done. Commented Apr 13, 2022 at 11:55
• In section 3.2 of that last reference by Six, a recursive procedure is given that looks like what I had imagined might exist. If it works it out I'll accept this as having provided a link to the answer. Commented Apr 13, 2022 at 17:40
• I implemented the recursive procedure, but was getting sign errors at intermediate steps that could be a bug, but haven't tracked it down. In the meantime, the solution below based on Macke et al. was provided and seems more straightforward, so I'm going with that. Commented Apr 19, 2022 at 15:58

This is not an additional answer but some summary from the results found in the references, particularly this one.

Let $${\bar {\bf z}} =a {\bar {\bf r}}+ b\, {\bf s} \, {\bar u} \tag 1$$

where $$a,b$$ are positive constants, $${\bar u}$$ is a vector of $$n$$ ones, $$s$$ and $${\bar {\bf r}}$$ are independent standard gaussian (scalar and $$n-$$multivariate resp). Then, $${\bar {\bf z}} \sim \mathcal N(0,\Sigma)$$ with $$\Sigma = a^2 I +b^2 U$$. Hence, by setting $$b=\sqrt{\tilde \rho}$$ and $$a = \sqrt{1-\tilde \rho}$$, $${\bar {\bf z}}$$ models our (positively) equicorrelated gaussian variables.

Let $${\bar {\bf x}}$$, where $${ \bf x_i}={\mathbf 1}_{z_i\ge 0}$$ be our correlated coins. And let $${ \bf y} =\sum_i { \bf x_i} \in 0,1\cdots n$$ be their sum.

Notice that $${\bar {\bf z}}$$ conditioned on $${\bf s}$$ are iid, with $$z_i | s \sim \mathcal N(bs, a^2)$$

Also $$x_i | s \sim \mathcal B_e(1 - \Phi(-sb/a))=B_e(\Phi(sb/a))$$ and $$y | s \sim \mathcal B(n,\Phi(sb/a))$$ , where $$\Phi$$ is the standard cumulative gaussian function, and $$\mathcal B_e$$, $$\mathcal B$$ denotes the Bernoulli and Binomial distributions.

Now, noting that $$H({ \bf y} | {\bar {\bf x}})=0$$, we get

\begin{align} H({\bar {\bf x}})&=H({ \bf y})+H({\bar {\bf x}}|{ \bf y})\\ &=H({ \bf y})+ \sum_y P({ \bf y}=y) H({\bar {\bf x}}|{ \bf y}=y) \\ &=H({ \bf y}) + \sum_y P({ \bf y}=y) \log \binom{n}{y}\\ &= \sum_y P({ \bf y}=y) \log \frac{\binom{n}{y}}{P({ \bf y}=y)} \tag 2\\ &= n-D( p({\bf y}) || {\mathcal B}(n,1/2)) \end{align}

And we can compute $$P({ \bf y}=y)$$ by integrating:

$$p({\bf y}) = \int_{-\infty}^{\infty} P(y|s) P(s) ds = \int_{-\infty}^{\infty} \mathcal B(n,\Phi(sb/a)) \, \phi(s) ds \tag 3$$

This already allows us to compute numerically the entropy. For large $$n$$, we can instead write

$$H({\bar {\bf x}}) = H({\bar {\bf x}}|s) + I({\bar {\bf x}}; s) \tag 4$$

where $$I()$$ is the mutual information (some care is needed here because $${\bar {\bf x}}$$ is discrete and $$s$$ is continuous, but the equation is justified nevertheless). The first term is linear in $$n$$ (conditioning on $$s$$ makes the components $$x_i$$ iid) and the second term is $$O(\log n)$$, because

$$I({\bar {\bf x}}; s) \le I({\bar {\bf z}}; s) = h(z)-h(z|s)= \frac12 (\log |\Sigma| - \log |a^2 I|) = \frac12 \log(1+b^2n) \tag 5$$

Then $$H({\bar {\bf x}}) = n \int \phi(s) h_b \left(\Phi(s \,b/a)\right) ds+ O(\log n) \tag 6$$

where $$h_b(p)=-p \log p -(1-p)\log(1-p)$$ is the binary entropy function.

The same approach can be used to map $$\tilde \rho \to \rho$$:

\begin{align} \rho &= 4 P(x_1=1,x_2=1)-1 \\ &= 4 \int P(x_1=1,x_2=1|s) p(s)ds -1 \\ &= 4 \int \Phi^2(sb/a) \phi(s) ds -1 \\ &= 4 \int \Phi^2\left(s \sqrt{\frac{\tilde \rho}{1-\tilde \rho}}\right) \phi(s)\, ds -1 \\ &= \frac{2}{\pi} \sin^{-1}({\tilde \rho}) \end{align}

• The last integral appears to be same performed in the OP and can be solved analytically using spherical coordinates giving $\tilde{rho}=\sin(\frac{\pi}{2}\rho)$. The claim here is that it doesn't not in fact depend on $n$, which seems correct. Commented Apr 19, 2022 at 15:18
• Can the mutual information approach be used to compute the large-n limit for the entropy of $y$? This entropy is $1/2\log(n)$ at $\tilde{\rho}=0$ (the large-$n$ limit of $B(n,1/2)$), has its max of $\log(n)$ at $\tilde{\rho}=1/2$, and is $1$ for $\tilde{\rho}=1$. In the approach here, the entropy of $y|s$ is now $O(\log(n))$, so the mutual information term is not negligible, but maybe still tractable. Or maybe there is another approach? Commented May 10, 2022 at 3:25
• Empirically, the large n result looks like it approaches the same form in $(H-1)/(\log(n+1)-1)$ as the result for $n=2$ (which in turn is trivially obtained by grouping the 01 and 10 sample in $H_{N=2}(\rho)$ to get $H_{N=2}(\rho)-(1-\rho)/2$. Thus the conjecture is that the limiting form is $(H_{N=2}(\rho)-(1-\rho)/2-1)/(\log(3)-1)(\log(n+1)-1)+1=(H_{N=2}(\rho)-(1-\rho)/2-1)/(\log(3)-1)\log(n)$ for large $n$. Commented May 10, 2022 at 13:32
• @puelmato Interesting , I don't have answers now - perhaps you can make that into a new question Commented May 10, 2022 at 16:43