On the Estimation of Discrete Distributions

1. Introduction

Let us consider a finite set of n situations that may occur according to a distribution ${\left(p_k\right)}_{k=0}^{n-1}$. The probability of each situation is generally determined by repeating the experiment independently and using Laplace’s operational definition to obtain the estimators

$${\widehat{p}}_k={\displaystyle \frac{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{o}\mathrm{c}\mathrm{c}\mathrm{u}\mathrm{r}\mathrm{r}\mathrm{e}\mathrm{n}\mathrm{c}\mathrm{e}\mathrm{s} \mathrm{o}\mathrm{f} \mathrm{s}\mathrm{i}\mathrm{t}\mathrm{u}\mathrm{a}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n} \mathrm{k}}{\mathrm{N}\mathrm{u}\mathrm{m}\mathrm{b}\mathrm{e}\mathrm{r} \mathrm{o}\mathrm{f} \mathrm{r}\mathrm{e}\mathrm{p}\mathrm{e}\mathrm{t}\mathrm{i}\mathrm{t}\mathrm{i}\mathrm{o}\mathrm{n}\mathrm{s} \mathrm{o}\mathrm{f} \mathrm{t}\mathrm{h}\mathrm{e} \mathrm{e}\mathrm{x}\mathrm{p}\mathrm{e}\mathrm{r}\mathrm{i}\mathrm{m}\mathrm{e}\mathrm{n}\mathrm{t}}}$$

In addition, the central limit theorem provides an estimate of the difference between the true value and the estimate [1].

This empirical method generally works quite well, but poses a few problems. On the one hand, this estimator and its uncertainty depend on the number of repetitions, giving good results only when this number is sufficiently large, which is a rather vague concept. On the other hand, it is difficult to generalize this definition when the observation of the situation is unclear. However, a situation is often determined using a measuring instrument that produces a noisy result.

The aim of this work is to base the estimation of probabilities associated with events on Bayes’ theorem, starting from Jeffreys’ non-informative priors.

2. Ingredients

2.1. Elements of Bayesian Theory of Statistical Decision-Making

This theory first distinguishes between two spaces: the space of states $\theta$ and the space of observations x[3]. Each observation is assumed to be made in the presence of an unknown state. The measurement model is the stochastic relation linking the state to the result of the measurement. For a given state, we give a transition probability1$\mathrm{p}(x\mid \theta )$, which gives the distribution of observations x given $\theta$. Given an a priori distribution $\mathrm{p}\left(\theta \right)$ over states and an observation x, Bayes’ theorem gives the a posteriori distribution over states

$$\mathrm{p}(\theta \mid x)={\displaystyle \frac{\mathrm{p}(x\mid \theta )\mathrm{p}\left(\theta \right)}{\int d\theta \mathrm{p}(x\mid \theta )\mathrm{p}\left(\theta \right)}}$$

When the state space is a differentiable manifold and the measurement model is twice differentiable with respect to $\theta$, we define a prior distribution called Jeffreys’ prior as follows [2]. First, we form the Fisher information matrix

$$I_{ij}\left(\theta \right)=-\int dx\mathrm{p}(x\mid \theta ){\displaystyle \frac{{\partial }^{2}ln\mathrm{p}(x\mid \theta )}{\partial {\theta }_i\partial {\theta }_j}}$$

Jeffreys’ prior is then defined as

$${\mathrm{p}}_J\left(\theta \right)d\theta =\sqrt{detI\left(\theta \right)}d\theta$$

A fundamental property of this distribution is its invariance. Indeed, it can be shown that ${\mathrm{p}}_J\left(\theta \right)d\theta$ does not depend on the choice of coordinate system.

2.2. Hyperspheres

The estimators constructed in this work make use of unit hyperspheres. Let us recall some of their properties [4].

In Euclidean space ${\mathbb{R}}^n$, the unit hypersphere of dimension $n-1$ is the differentiable manifold

$${\mathbb{S}}^{n-1}=\{\psi \in {\mathbb{R}}^n\mid \sum _k{\psi }_k^2=1\}$$

The use of the letter $\psi$ will be justified below. The usual spherical coordinates are given by

$$\begin{gathered} {\psi }_0=sin{\theta }_{n-2}sin{\theta }_{n-3}\dots sin{\theta }_{1}sin{\theta }_0 \\ {\psi }_1=sin{\theta }_{n-2}sin{\theta }_{n-3}\dots sin{\theta }_{1}cos{\theta }_0 \\ {\psi }_2=sin{\theta }_{n-2}sin{\theta }_{n-3}\dots cos{\theta }_1 \\ \dots \\ {\psi }_{n-1}=cos{\theta }_{n-2} \end{gathered}$$

(1)

and the normalized invariant measure under $SO\left(n\right)$ is given by

$$d\mu \left(\psi \right)={\displaystyle \frac{1}{A_{n-1}}}{sin}^{n-2}{\theta }_{n-2}{sin}^{n-3}{\theta }_{n-3}\dots sin{\theta }_{1}d{\theta }_{0}d{\theta }_1\dots d{\theta }_{n-2}$$

(2)

$$A_{n-1}={\displaystyle \frac{2{\pi }^{n/2}}{\Gamma (n/2)}}$$

(3)

$$\Gamma (1/2)=\sqrt{\pi }$$

(4)

Let $\mathit{m}=(m_0\dots m_{n-1})$ be a sequence of positive integers. The even moments relative to this measure are given by

$$H\left(\mathit{m}\right)=\int {\psi }_0^{2m_0}\dots {\psi }_{n-1}^{2m_{n-1}}d\mu \left(\psi \right)$$

(5)

We can show that the function H above is a generalization of the beta function

$$H\left(\mathit{m}\right)={\displaystyle \frac{\Gamma (n/2)}{{\pi }^{n/2}}}{\displaystyle \frac{{\prod }_k\Gamma (m_k+1/2)}{\Gamma \left({\sum }_k(m_k+1/2)\right)}}$$

(6)

The normalization factor is chosen so that $H\left(0\right)=1$

3. Probability Estimators for Known Situations

The problem is estimating a distribution $p_k$. This distribution can be viewed as a point in the $n-1$ dimensional simplex

$${\Delta }^{n-1}=\{\mathbf{p}\mid 0⩽p_k⩽1,\sum _k{p}_k=1\}$$

This simplex will be the state space for a measurement model whose observations are situations. In a canonical way, the probability of obtaining situation k given that the distribution is $\mathbf{p}$ is obviously

$$\mathrm{p}(k\mid \mathbf{p})=p_k$$

This canonical measurement model is twice differentiable with respect to $\mathbf{p}$. We can therefore calculate the associated Jeffreys prior ${\mathrm{p}}_J\left(\mathbf{p}\right)$. It was stated above that this distribution is independent of the parameterization of the simplex. Let us choose the new coordinates

$${\psi }_k=\sqrt{p_k}$$

(7)

Since ${\sum }_k{\psi }_k^2={\sum }_k{p}_k=1$, we obtain a measure model on a state space that is the positive part of the hypersphere ${\mathbb{S}}^{n-1}$

$$\mathrm{p}(k\mid \psi )={\psi }_k^2$$

Now, it turns out that Jeffreys’ prior for this model is the invariant measure $\mu$ mentioned above:

$$d{\mathrm{p}}_J\left(\psi \right)=d\mu \left(\psi \right)$$

(8)

This result allows us to apply Bayes’ theorem to a sequence of independent observations. Let

$$\mathit{i}=(i_0\dots i_{t-1})$$

be such a sequence of length t, we obtain the posterior on the sphere ${\mathbb{S}}^{n-1}$ given $\mathit{i}$

$$d\mathrm{p}(\psi \mid \mathit{i})={\displaystyle \frac{{\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}{\int {\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}}$$

(9)

Introducing the function $\mathit{m}$ giving the multiplicity of an index k in the sequence $\mathit{i}$:

$$\mathit{m}{\left(\mathit{i}\right)}_k=Card\{s\mid i_s=k\}$$

we obtain (implying the sequence $\mathit{i}$ to simplify the notation)

$$d\mathrm{p}(\psi \mid \mathit{i})={\displaystyle \frac{{\prod }_k{\psi }_k^{2m_k}d\mu \left(\psi \right)}{H\left(\mathit{m}\right(\mathit{i}\left)\right)}}$$

(10)

The probability of occurrence of a situation k can be estimated by the posterior mean,

$${\widehat{p}}_k=\mathbb{E}\left(p_k\right)={\displaystyle \frac{\int {\psi }_k^2{\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}{H\left(\mathit{m}\right(\mathit{i}\left)\right)}}$$

(11)

The associated uncertainties will be given by the covariance matrix

$$\widehat{p_k{p}_l}=\mathbb{E}\left(p_k{p}_l\right)={\displaystyle \frac{\int {\psi }_k^2{\psi }_l^2{\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}{H\left(\mathit{m}\right(\mathit{i}\left)\right)}}$$

(12)

Let us define ${\mathit{e}}_k$ as a sequence of integers that is zero except for the $k^{th}$, whose value is 1. It is easy to see that

$${\widehat{p}}_k={\displaystyle \frac{H(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k)}{H\left(\mathit{m}\right(\mathit{i}\left)\right)}}$$

(13)

$$\widehat{p_k{p}_k}={\displaystyle \frac{H(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k+{\mathit{e}}_l)}{H\left(\mathit{m}\right(\mathit{i}\left)\right)}}$$

(14)

The fact that ${\sum }_h{m}_h=t$ and the properties of the Gamma function [5] immediately give

$${\widehat{p}}_k={\displaystyle \frac{m_k+1/2}{t+n/2}}$$

(15)

$$\widehat{p_k^2}={\displaystyle \frac{(m_k+1/2)(m_k+3/2)}{(t+n/2)(t+n/2+1)}}$$

(16)

$$\widehat{p_k{p}_l}={\displaystyle \frac{(m_k+1/2)(m_l+1/2)}{(t+n/2)(t+n/2+1)}}$$

(17)

4. Self-Consistent Random Sequences

Let $\mathit{i}$ be a sequence of length $t-1$. Equation (15) therefore gives the Bayesian estimates for the $p_k$ if we have observed the sequence $\mathit{i}$, which we will also denote by

$${\widehat{p}}_k={\mathrm{p}}_J(k\mid i_{t-2},i_{t-3}\dots i_0)$$

At t, let us randomly draw a new situation $i_{t-1}$ according to these $p_k$. Let us form a new sequence by concatenation: $(\mathit{i},i_{t-1})$. Nothing prevents us from repeating the process for this new sequence.

We say that a sequence $\mathit{i}$ is constructed in a self-consistent manner if

• $i_0$ is drawn uniformly from the n situations
• at each step, $i_{s-1}$ is drawn according to the probability $p_J(k\mid i_{s-2}\dots i_0)$

We note that the probability of obtaining a sequence $\mathit{i}=(i_0\dots i_{t-1})$ in this way is

$$p_{ac}\left(\mathit{i}\right)=p_J(i_{t-1}\mid i_{t-2}\dots i_0)p_J(i_{t-2}\mid i_{t-3}\dots i_0)\dots p_J(i_1\mid i_0)p_J\left(i_0\right)$$

(18)

Now, it turns out that this probability is given by the function $H\circ \mathit{m}$

$$p_{ac}\left(\mathit{i}\right)=H\left(\mathit{m}\left(\mathit{i}\right)\right)$$

(19)

Since no additional assumptions are made, we will call the probability $p_{ac}$a priori (non-informative) self-consistent on sequences $\mathit{i}$ of given length t. Sampling according to this prior will be useful later for estimation in uncertain situations.

Note

This way of generating a random sequence of numbers provides a possible answer to the question: Given a sequence, what is its probability of occurrence? It can be easily seen that the sequences with maximum self-consistent probability are also those with minimum entropy.

Incidentally, for sequences of given length t, we have

$$\sum _{\mathit{i}}H\left(\mathit{m}\left(\mathit{i}\right)\right)=1$$

From relations (18) and (19), we also have

$$H(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k)=p_J(k\mid \mathit{i})H\left(\mathit{m}\left(\mathit{i}\right)\right)H$$

(20)

$$(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k+{\mathit{e}}_l)=p_J(k\mid \mathit{i})p_J(l\mid (\mathit{i},k))H\left(\mathit{m}\left(\mathit{i}\right)\right)$$

(21)

5. Generalization to Fuzzy Situations

In reality, situations are often measured by noisy instruments. Instead of knowing the situation at a stage s, the observer only has a measurement result x from a measurement model $\mathrm{p}(x\mid k)$. For a given x, this vector, whose indices are the situations, is often called the likelihood function. For a sequence of steps $s=0\dots t-1$, the observations2$(x_0\dots x_{t-1})$ give a likelihood matrix whose time is the row index and the situation is the column index

$${\ell }^s\left(k\right)=\mathrm{p}(x_s\mid k)$$

(22)

Let’s return to the situation on the hypersphere. Suppose that at time s, the prior probability on ${\mathbb{S}}^{n-1}$ is $d{\eta }_s\left(\psi \right)$. The appearance of a result $x_s$ for the measurement model $\mathrm{p}(x\mid k)$, combined with the canonical model $\mathrm{p}(k\mid \psi )={\psi }_k^2$, gives, by Bayes’ theorem, an a posteriori that will be taken as a priori in the next step:

$$d{\eta }_{s+1}\left(\psi \right)={\displaystyle \frac{{\sum }_k{\ell }^s\left(k\right){\psi }_k^{2}d{\eta }_s\left(\psi \right)}{{\sum }_k{\ell }^s\left(k\right)\int {\psi }_k^{2}d{\eta }_s\left(\psi \right)}}$$

(23)

We can therefore write

$$d{\eta }_t\left(\psi \right)={\displaystyle \frac{{\sum }_{\mathit{i}}{\ell }^0\left(i_0\right)\dots {\ell }^{t-1}\left(i_{t-1}\right){\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}{{\sum }_{\mathit{i}}{\ell }^0\left(i_0\right)\dots {\ell }^{t-1}\left(i_{t-1}\right)\int {\psi }_{i_0}^2\dots {\psi }_{i_{t-1}}^{2}d\mu \left(\psi \right)}}$$

(24)

Let

$$\ell \left(\mathit{i}\right)=\prod _s{\ell }^s\left(i_s\right)$$

and use the definition of the function H to obtain the estimators at $t-1$

$$\mathbb{E}\left(p_k\right)={\displaystyle \frac{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k)}{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}}$$

(25)

$$\mathbb{E}\left(p_k{p}_l\right)={\displaystyle \frac{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H(\mathit{m}\left(\mathit{i}\right)+{\mathit{e}}_k+{\mathit{e}}_l)}{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}}$$

(26)

which can also be written as

$$\mathbb{E}\left(p_k\right)={\displaystyle \frac{{\sum }_{\mathit{i}}{p}_J(k\mid \mathit{i})\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}}$$

(27)

$$\mathbb{E}\left(p_k{p}_l\right)={\displaystyle \frac{{\sum }_{\mathit{i}}{p}_J(k\mid \mathit{i})p_J(l\mid (\mathit{i},k))\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}{{\sum }_{\mathit{i}}\ell \left(\mathit{i}\right)H\left(\mathit{m}\left(\mathit{i}\right)\right)}}$$

(28)

The preceding equations suggest a new use for Bayes’ theorem. Indeed, given the self-consistent prior $p_{ac}$ on sequences of length t and the product measure model giving the likelihood function ℓ, we can form the posterior on the sequences

$$p_{ac}(\mathit{i}\mid \ell )={\displaystyle \frac{\ell \left(\mathit{i}\right)p_{ac}\left(\mathit{i}\right)}{{\sum }_{\mathit{j}}\ell \left(\mathit{j}\right)p_{ac}\left(\mathit{j}\right)}}$$

(29)

and we obtain

$$\mathbb{E}\left(p_k\right)=\sum _{\mathit{i}}{p}_J(k\mid \mathit{i})p_{ac}(\mathit{i}\mid \ell )$$

(30)

$$\mathbb{E}\left(p_k{p}_l\right)=\sum _{\mathit{i}}{p}_J(k\mid \mathit{i})p_J(l\mid (\mathit{i},k))p_{ac}(\mathit{i}\mid \ell )$$

(31)

Equations (30) and (31) are integrals that can be calculated using Monte Carlo. It is possible to sample the posterior $p_{ac}(\mathit{i}\mid \ell )$ as follows:

• draw $i_0$ according to the weights ${\ell }^0\left(k\right)$
• draw $i_1$ according to the weights ${\ell }^1\left(k\right)p_J(k\mid i_0)$
• …
• draw $i_{t-1}$ according to the weights ${\ell }^{t-1}\left(k\right)p_J(k\mid i_{t-2}\dots i_0)$

A sample of r sequences ${\mathit{i}}_s$ of this type gives the MC estimators

$$\widehat{\mathbb{E}\left(p_k\right)}={\displaystyle \frac{1}{r}}\sum _s{p}_J(k\mid {\mathit{i}}_s)$$

(32)

$$\widehat{\mathbb{E}\left(p_k{p}_l\right)}={\displaystyle \frac{1}{r}}\sum _s{p}_J(k\mid {\mathit{i}}_s)p_J(l\mid ({\mathit{i}}_s,k))$$

(33)

as well as MC uncertainties

$$u_{\mathbb{E}\left(p_k\right)}={\displaystyle \frac{1}{\sqrt{r}}}\sqrt{{\displaystyle \frac{1}{r}}{\sum }_s{p}_J{(k\mid {\mathit{i}}_s)}^2-{\widehat{\mathbb{E}\left(p_k\right)}}^2}$$

(34)

$$u_{\mathbb{E}\left(p_k{p}_l\right)}={\displaystyle \frac{1}{\sqrt{r}}}\sqrt{{\displaystyle \frac{1}{r}}{\sum }_s{p}_J{(k\mid {\mathit{i}}_s)}^2{p}_J{(l\mid ({\mathit{i}}_s,l))}^2-{\widehat{\mathbb{E}\left(p_k{p}_l\right)}}^2}$$

(35)

Note

Simulating this sampling is very easy, for example in C++, and gives excellent results for likelihood matrices with 3 situations and 10,000 observation times with an r of around one million. The advantage of sampling on the posterior is that it directly locates the sequences that contribute significantly to the overall sum. Recall that the sums (30) and (31) contain $n^t$ terms.

6. Conclusions

Replacing $p_k$ with ${\psi }_k=\sqrt{p_k}$ makes Jeffreys’ prior uniform. In general, probabilistic quantities can be parameterized by the square root, which is their natural expression. This phenomenon is found in quantum mechanics [6], where probability is expressed as the square of the modulus of a wave function, hence the choice of the Greek letter $\psi$ for the points of the hypersphere.

The generalization of beta functions given by the function H bridges the gap between continuous calculus on wave functions and sums over sequences, which are computable by MC.

An application of this work is under development for stationary Markov processes.

7. Summary

7.1. Seen by $\psi \in {\mathbb{S}}^{n-1}$

$$d{\eta }^0\left(\psi \right)=d\mu \left(\psi \right)\mathrm{i}\mathrm{n}\mathrm{v}\mathrm{a}\mathrm{r}\mathrm{i}\mathrm{a}\mathrm{n}\mathrm{t}$$

(36)

$$d{\eta }^t\left(\psi \right)=d{\eta }^{t-1}(\psi \mid {\ell }^{t-1})={\displaystyle \frac{{\sum }_k{\ell }^{t-1}\left(k\right){\psi }_k^{2}d{\eta }^{t-1}\left(\psi \right)}{{\sum }_k{\ell }^{t-1}\left(k\right)\int {\psi }_k^{2}d{\eta }^{t-1}\left(\psi \right)}}$$

(37)

$${\mathbb{E}}^t(p_k\mid \ell )=\int {\psi }_k^{2}d{\eta }^t\left(\psi \right)$$

(38)

$${\mathbb{E}}^t(p_k{p}_l\mid \ell )=\int {\psi }_k^2{\psi }_l^{2}d{\eta }^t\left(\psi \right)$$

(39)

7.2. Seen by $\mathit{i}\in {\{0\dots n-1\}}^{\mathbb{N}}$

$${\mathit{i}}^t=(i_0\dots i_{t-1})$$

(40)

$$p_{ac}^t\left(\mathit{i}\right)=p_J(i_{t-1}\mid {\mathit{i}}^{t-1})p_{ac}^{t-1}$$

(41)

$$\left(\mathit{i}\right){\ell }^t\left(\mathit{i}\right)=\prod _{s=0}^{t-1}{\ell }^s\left(i_s\right)$$

(42)

$$p_{ac}^t(\mathit{i}\mid \ell )={\displaystyle \frac{{\ell }^t\left(\mathit{i}\right)p_{ac}^t\left(\mathit{i}\right)}{{\sum }_{{\mathit{i}}^t}{\ell }^t\left(\mathit{i}\right)p_{ac}^t\left(\mathit{i}\right)}}$$

(43)

$${\mathbb{E}}^t(p_k\mid \ell )=\sum _{{\mathit{i}}^t}{p}_J(k\mid {\mathit{i}}^t)p_{ac}^t(\mathit{i}\mid \ell )$$

(44)

$${\mathbb{E}}^t(p_k{p}_l\mid \ell )=\sum _{{\mathit{i}}^t}{p}_J(l\mid {\mathit{i}}^{t}k)p_J(k\mid {\mathit{i}}^t)p_{ac}^t(\mathit{i}\mid \ell )$$

(45)