Bayes estimator for multinomial parameters and Bhattacharyya distances
Christopher Ferrie
Centre for Engineered Quantum Systems, School of Physics, The University of Sydney, Sydney, NSW, Australia

Robin Blume-Kohout

arXiv:1612.07946v1 [math.ST] 23 Dec 2016

Sandia National Laboratories, Albuquerque, New Mexico, 87185
(Dated: December 26, 2016)
We derive the Bayes estimator for the parameters of a multinomial distribution under two loss
√
functions (1 − B and 1 − B2 ) that are based on the Bhattacharyya coefficient B( p, q) = ∑ pk qk .
We formulate a non-commutative generalization relevant to quantum probability theory as an open
problem. As an example application, we use our solution to find minimax estimators for a binomial
parameter under Bhattacharyya loss (1 − B2 ).

I.

INTRODUCTION

In statistical decision theory, a Bayes estimator for a prior µ is a point estimator whose expected risk (with respect
to µ) is as small as possible [1]. Bayes estimators have utility beyond strictly Bayesian statistics; they also provide
lower bounds on minimax risk (a frequentist concept), and the greatest lower bound (That is, the maximum Bayes
risk over all priors) coincides with the minimax risk. This duality enables efficient numerical algorithms to find
minimax estimators by optimizing over Bayesian priors, rather than over the computationally intractable space of
estimators [2].
Different loss functions produce different risk functions, and the Bayes estimator is sensitive to this choice. For a
certain class of risks known as Bregman divergences the Bayes estimator is always the mean of the posterior distribution [3]. For the multinomial model that we consider, both mean-squared-error and Kullback-Leibler divergence
(relative entropy) are Bregman divergences. But not every useful risk function has this property. We consider two
loss functions defined through the Bhattacharyya coefficient B( p, q), a Fisher-adjusted measure of distinguishability
between probability distributions that finds use in machine learning (as a metric on probability distributions [4]) and
in quantum information theory (as fidelity [5, 6], the most commonly used measure of quantum state distinguishability in both theoretical and experimental studies).
In this paper, we find the Bayes estimators for 1 − B and 1 − B2 by mapping their respective optimizations to
simple linear programs. We also state, but do not solve, the non-commutative generalization of the problem, which
has applications in quantum information theory. Finally, we demonstrate the utility of the Bayes estimator by finding
the minimax estimator for a binomial parameter.
II.

PROBLEM SET-UP, REVIEW AND CHOICE OF LOSS FUNCTIONS

Suppose we have a coin, a K-sided die, or any process that generates i.i.d. samples from a finite set of size K. We
wish to estimate the probabilities { p0 , . . . pK −1 }. They form a parameter vector p which belongs to the K-simplex
(
)

4K =

( p0 , p1 . . . , pK−1 ) pk ≥ 0 ∀k, ∑ pk = 1 .

(1)

k

Estimation requires some data. From the data, we wish to produce an estimate of p, call it p̂. Although the data
are most commonly obtained by repeated sampling directly from p, other experiments are possible. We might have
“noisy” data (drawn from an ancillary distribution q = Σp, where the stochastic matrix Σ is known [7]), or side
information. Happily, the Bayes analysis turns out to be independent of how the data were obtained. Instead, it
depends only on the posterior probability density Pr( p|data)—which is to say that, given the posterior, the expected
risk is independent of the likelihood Pr(data| p).
The quality of an estimate p̂ is formalized by a loss function L( p, p̂) that quantifies how bad it is to report an
estimate p̂ when the truth is p. The best estimate should be p itself, so L( p, p) ≤ L( p, p̂) for all p̂. A “good”
estimator is one that minimizes expected loss, or risk:
r p = Edata [ L( p, p̂(data)].

(2)

2
The risk is a function of the underlying p, and it is generally impossible to minimize risk for all p simultaneously.
The Bayesian solution is to consider (and minimize) the average risk with respect to a given prior µ( p). This is called
the Bayes risk of the estimator p̂ with respect to µ:
r ( p̂) = E p,data [ L( p, p̂(data)].

(3)

An estimator p̂B that minimizes the Bayes risk is called a Bayes estimator for µ( p). The Bayes estimator can also be
defined (under normal circumstances; [8]) explicitly for each dataset as the minimizer of posterior risk. This identity
is extremely useful, because it means that a Bayes estimator can be constructed explicitly for any given value of data
without considering others.
Different loss functions yield different Bayes estimators. For example, for quadratic loss functions of the form
L( p, p̂) = ( p − p̂)T A( p − p̂),

(4)

where A > 0 is a positive definite matrix, the Bayes estimator is the mean of the posterior distribution [8]
p̂B (data) = E p|data [ p].

(5)

Another important loss function is the relative entropy or Kullback-Leibler divergence:
L( p, p̂) =

p

∑ pk log p̂kk .

(6)

k

Here, again, the Bayes estimator is the posterior mean [9]. This is not an accident; it holds true whenever the loss is
a Bregman divergence [3]. Quadratic loss and Kullback-Leibler loss are both Bregman divergences.
Unfortunately, not every useful and important loss function is a Bregman divergence. The loss functions we
consider here are not Bregman divergences, and their Bayes estimators are not the posterior mean. They are based
on the Bhattacharyya coefficient [10],
p
B( p, p̂) = ∑ pk p̂k ,
(7)
k

which quantifies the statistical indistinguishability or fidelity of two distributions. In different contexts, both L1 =
1 − B and L2 = 1 − B2 are useful loss functions, and so we find Bayes estimators for both of them. In the limit of
small deviations ( p̂ is close to p), they coincide up to a factor of 2. Both are Fisher-adjusted, meaning that their 2nd
order series expansions are proportional to the Fisher metric. However, L1 and L2 have different global behavior,
and therefore different Bayes estimators.
Since the minimizer of 1 − B p is the maximizer of B p , we define two Bayes estimators as follows:
p̂1 = argmax E p [ B( p, p̂)],

(8)

p̂2 = argmax E p [ B2 ( p, p̂)],

(9)

p̂
p̂

subject to p̂ ∈ 4K . Note that, since the Bayes estimator depends only on the posterior, with no explicit dependence
on the prior or the data, we hereafter drop the conditional notation on the data. That is, E p means expectation with
respect to an arbitrary distribution of p.
III.

BAYES ESTIMATORS

Let us now derive expressions for the Bayes estimators defined in Equations 8 and 9. A very useful way of writing
√
B is to define p as the element-wise square root of p, then write:
B( p, p̂) =

√

p·

p

p̂.

The posterior risk involves an average of B over p. Both this average and the dot product in B are linear in
means that they commute, so:
p
√
E p [ B( p, p̂)] = E p [ p] · p̂.

(10)

√

p, which
(11)

3
This also implies the normalization
1=

√ √
√
p k p k = k p k2 .

∑ pk = ∑
k

(12)

k

Thus, the problem in Equation 8 becomes
p
√
maximize E p [ p] · p̂
p p
subject to
p̂ · p̂ = 1,
p
p̂ ≥ 0.
This is a textbook [11] convex optimization problem, and the solution is
√
p
E p [ p]
,
p̂1 =
√
kE p [ p]k

(13)

(14)

implying

√
E p [ p ]2
p̂1 =
√
√ ,
(E p [ p]) · (E p [ p])

(15)

where the square in the numerator is element-wise.
Minimizing 1 − B2 , as in Equation 9, requires a bit more work, but we can again reduce it to a linear program. We
start by defining
P=
the projection operator onto

√

√ √ T
p p ,

(16)

p, and similarly for P̂. The normalization condition on P is
Tr[ P] = 1.

(17)

In terms of these objects, the squared Bhattacharyya coefficient is
B2 ( p, p̂) = Tr[ P P̂].

(18)

This expression is now linear in P, which allows us to apply the trick of commuting the expectation through:
E p [ B2 ( p, p̂)] = Tr[E p [ P] P̂].
The matrix elements of the expectation on the RHS can be written out explicitly as

p
E p [ P] ij = E p [ pi p j ].

(19)

(20)

Putting everything together, Equation 9 becomes
maximize

Tr [E p [ P] P̂]

subject to

Tr[ P̂] = 1,

(21)

P̂ ≥ 0,
where P̂ ≥ 0 means positive semi-definiteness as a matrix. This is another textbook [11] linear semidefinite program.
If we define a as the normalized eigenvector of E p [ P] with maximal eigenvalue, then the solution to Equation 21 is
P̂2 = aaT .

(22)

Explicitly, as a vector in 4K , we have
p̂2 = a2 ,
(23)
√
where a = ( a0 , a1 , . . . , aK −1 ) is the eigenvector of E p [ pi p j ] with maximal eigenvalue, and the square is elementwise as in Eq. 15.

4
IV.

NON-COMMUTATIVE GENERALIZATION: A QUANTUM (OPEN) PROBLEM

Quantum mechanics can be described as a non-commutative generalization of probability theory [12]. Quantum
random variables (That is, the “state” of a quantum system) are represented not by K-element vectors of probabilities
p, but by K × K self-adjoint complex matrices or density matrices. In this framework, observable events are also
represented by K × K self-adjoint matrices called effects. States can be sampled or observed in more than one way.
The quantum analogue to “sampling” is “measuring” the system, and the observer must choose how to measure
it. This means specifying a mutually exclusive and exhaustive set of possible events, each represented by a positive
semidefinite d × d matrix. The set of event matrices { Ek } is called a positive operator valued measure (POVM), and
satisfies the condition ∑k Ek = 1l. If an quantum system is state ρ is observed, then the probability of observing
the event represented by effect E is given by Pr( E|ρ) = Tr(ρE). As a result, the law of total probability (“exactly
one of the possible events occurs”) manifests itself as the constraint Trρ = 1, and the non-negativity of probabilities
corresponds to ρ ≥ 0. Thus, the set of quantum states (analogous to the simplex ∆ of valid probability vectors) is

S = {ρ|ρ ≥ 0, Tr(ρ) = 1} .

(24)

The subset of S containing only diagonal matrices (“classical” states) is in 1:1 correspondence with K-element probability vectors, since ρ ≥ 0 ensures positivity of the diagonal entries, and Tr(ρ) = 1 ensures their normalization.
Thus, if we consider only states that are simultaneously diagonalizable, we can identify p ≡ ρ and p̂ ≡ ρ̂, and the
Bhattacharyya coefficient has a simple expression in terms of ρ and ρ̂:
p 
ρρ̂ .
(25)
B( p, p̂) = Tr
Its square—the B2 risk that we considered previously—is called classical fidelity in the quantum mechanics literature.
Hence, our solution to Equation 21 can be applied directly to estimation of quantum states, whenever the support of
the averaging measure (e.g. a posterior distribution) consists entirely of commuting quantum states:
p
maximize Eρ [Tr( ρρ̂)2 ]
(26)
subject to Tr(ρ̂) = 1,
ρ̂ ≥ 0.
However, quantum states are not usually restricted to be diagonal, which means that in the estimation of quantum
states (quantum tomography [13]), the posterior will have support on non-commuting density matrices. This requires
new definitions of loss and risk. Without getting into the details (see Ref. [6]), the general definition of quantum
fidelity between two states is:

 q
√ √ 2
F (ρ, ρ̂) = Tr
ρρ̂ ρ ,
(27)
and the corresponding loss function is infidelity or L = 1 − F. For co-diagonal states ρ̂ and ρ, fidelity coincides with
B2 . Thus, if the posterior does happen to be supported only on mutually commuting states, then our solution to
Equation 21 gives the Bayes estimator. In general, however, the problem of finding a Bayes estimator for quantum
infidelity is:
maximize Eρ [ F (ρ, ρ̂)]
subject to

Tr(ρ̂) = 1,
ρ̂ ≥ 0.

(28)

This problem appears to be quite difficult (some bounds are given in [14]), and we do not attempt to solve it in full
generality here.
V.

EXAMPLE: BINOMIAL PARAMETER
A.

Canonical beta prior

As an interesting and useful application of our main result, we now compute (and examine) the Bayes estimator
p̂2 for L2 = 1 − B2 in the simple case of a binomial distribution—that is, a directly observed coin toss or Bernoulli

5
process. Here, p = ( p0 , p1 ), with p1 = 1 − p0 . The canonical (conjugate) prior is a beta distribution:
β −1 β −1

Pr( p) =

p0 p1
dp,
Beta( β, β)

(29)

where Beta() is the beta function. If the coin is flipped N times, yielding n heads and N − n tails, then the posterior
is
n + β −1 N − n + β −1

Pr( p|n) =

p0
p1
dp.
Beta(n + β, N − n + β)

(30)

√
The matrix E p [ p j pk ] is given explicitly in terms of beta functions by


Beta( β + n + 1, β + N − n)
Beta( β + n + 1/2, β + N − n + 1/2)
Beta( β + n + 1/2, β + N − n + 1/2)
Beta( β + n, β + N − n + 1)


/Beta( β + n, β + N − n).

(31)

To compute the Bayes estimator p̂2 , we diagonalize this matrix, extract the eigenvector corresponding to the larger
eigenvalue, and square its entries. This can be written down in closed form, but the result is lengthy, messy, and
opaque. Instead, we illustrate its properties for N = 10 flips and for conjugate priors with β = 1 (Laplace’s prior)
and β = 1/2 (Jeffrey’s prior), in Fig. 1. For comparison, we also show (1) the posterior mean estimator, and (2) the
maximum likelihood estimator.

Jeffreys Prior (β = 1/2)

0.5

Laplace Prior (β = 1)

0.5

0.3

0.3

p̂

0.4

p̂

0.4

0.2

0.2

MLE
Mean
Bayes

0.1
0.0

0

1

2

n

3

4

MLE
Mean
Bayes

0.1
5

0.0

0

1

2

n

3

4

5

FIG. 1: The maximum likelihood estimator (MLE), optimal Bayes estimator and the posterior mean for N = 10 coin flips. On
the left, the latter two estimators begin with a canonical beta prior with β = 1/2 (Jeffrey’s prior). On the right, the same but with
β = 1 (Laplace’s uniform prior). The fidelity performance of these estimators is shown in Figure 3.

The Bayes estimator is generally similar to the posterior mean estimator. Both “hedge” away from p = 0 and
p = 1, in contrast to the MLE, but the Bayes estimator is more aggressive (less hedged) than the mean. Since the
mean estimator is Bayes for Bregman divergences, this implies that the 1 − B2 loss function is more forgiving of
extreme estimates (That is, ones close to p = 0 or p = 1) than any Bregman divergence.
However, although the difference between the two estimators (Bayes and posterior mean) appears significant, the
difference in their performance is not. While the Bayes estimator does have lower Bayes risk, the difference is small
even for N = 10. Figure 2 shows the relative suboptimality of the mean estimator, defined as:

R=

rBayes − rmean
,
rBayes

(32)

where rmean = E p|n [1 − B2 ( p, p̂mean (n))] and likewise for rBayes . The relative suboptimality of the mean estimator—
even at N = 10—is always less than 10−3 , decreases rapidly as a function of N and is largely independent of β.
Given that the mean estimator is far easier to calculate (the Bayes estimator requires computing the matrix
√
E p [ p j pk ] and its eigenvectors), we suggest that in practice for most purposes, the mean estimator is a completely
satisfactory heuristic. For small N, there is a noticeable difference in the estimate itself (see Fig. 3), but not in performance (see Fig. 2).

6

0.025

Relative Loss

0.020
0.015
0.010
0.005

1.1

1.0

0.9

0.8

β

0.7

0.6

0.5 100

90

80

70

60

50

40

30

20

0.000
10

N

FIG. 2: The relative difference, R in Eq. (32), between the average loss of the Bayes estimator rBayes and that of the mean estimator
rmean .

B.

Least favorable prior

Frequentist analyses of estimators average the risk only over the data (That is with respect to Pr(data| p)), rather
than over the joint distribution Pr(data, p). This more restricted average defines pointwise risk,
R( p, p̂) = Edata| p [ L( p, p̂)],

(33)

which retains a dependence on p that has to be gotten rid of somehow in order to define an “optimal” estimator.
Instead of averaging over a prior (to get the Bayes risk in Eq. (3)), the most common frequentist way to remove this
dependence is to consider the “worst case”, and maximize the risk over p:
Rmax ( p̂) = max R( p, p̂).
p

(34)

An estimator that minimizes the maximum risk is called a minimax estimator:
p̂minimax = arg min Rmax ( p̂).

(35)

p̂

The minimax criterion seeks an estimator whose performance is “pretty good” for any value of the parameters, without reference to prior probability. In practice, minimax estimators achieve (approximately) equal risk for all p. The
Bayes estimators for the binomial parameter that we derived in the previous section are not minimax; Figure 3 shows
that their pointwise risk varies with p0 . However, they’re pretty close—the variations are small (e.g., compared with
those of the MLE).
However, while the minimax estimator is unique, each prior has a different Bayes estimator. We can rewrite Eq. (3),
making this explicit, letting µ denote the prior density:
r (µ, p̂) = E p [ R( p, p̂)]

(36)

Somewhat remarkably, there is (at least under weak regularity conditions [1]) a prior whose Bayes estimator is
minimax. The least favorable prior is the prior whose Bayes estimator has the highest Bayes risk (out of all priors), and
its Bayes estimator is minimax. Explicitly, the estimator is
p̂LFP = arg max min r (µ, p̂)
µ

p̂

(37)

This remarkable relation between Bayesian and frequentist optimality is called Bayes-minimax duality:
p̂LFP = p̂minimax .

(38)

7

Jeffreys Prior (β = 1/2)

MLE
Mean
Bayes

Average 1 − B 2

0.04
0.03
0.02
0.01
0.00
0.0

0.1

0.2

p

0.3

0.4

0.09
0.08
0.07
0.06
0.05
0.04
0.03
0.02
0.01
0.00
0.0

Laplace Prior (β = 1)

MLE
Mean
Bayes

Average 1 − B 2

0.05

0.5

0.1

0.2

p

0.3

0.4

0.5

FIG. 3: The pointwise risk (1 − B2 ) of the estimators shown in Figure 1, as a function of p0 .

One useful consequence of Bayes-minimax duality is that it provides a method for finding minimax estimators by
searching for the least favorable prior. While the space of priors is large, it is enormously smaller than the space of
estimators. However, this approach is only practical if we have a closed-form expression for Bayes estimators, since
this allows fast and efficient computation of the Bayes risk for each prior considered. Our result enables this sort of
analysis for 1 − B and 1 − B2 loss functions.
We used our result to construct minimax estimators for the binomial parameter. Unlike the Bayes estimator,
the minimax estimators do not have explicit closed forms; the least favorable priors are not elegant, and require
numerical optimization. We performed two numerical optimizations. First, we did a restricted optimization over
conjugate priors (beta distributions, of the form given in Eq. 29), to find the minimax value of β. Then, we performed
an unrestricted optimization over all priors to find the true minimax estimator
In the first (restricted) optimization over conjugate priors, we found that the optimal β varied only weakly with
N, and was given by β ≈ 0.44 for all N. This is extremely consistent with other answers to the general question
“What value of β works best?”; a variety of other analyses have yielded answers close to the β = 1/2 that defines
Jeffrey’s estimator. Just for reference, we considered what prior would yield the lowest maximum risk if we were to
use the posterior mean estimator rather than the Bayes estimator derived in this paper. For this ad-hoc procedure,
the optimal β is approximately β ≈ 0.26. This estimator is not minimax in any sense, but it does support our general
observation that the Bayes estimator can be approximated by the posterior mean without much damage (0.26 and
0.44 are not very different, and the achieved maximum risks are also not very different).
To find the true minimax estimator, we used the algorithm of Kempthorne [2] to find the least favorable prior.
Roughly speaking, this algorithm iteratively adds support points to a discrete prior, each time maximizing the Bayes
risk over the small finite search space of support points and weights. It proceeds until the average risk and the
maximum risk are within a pre-defined tolerance. We provide pseudo-code for our implementation of the algorithm
in Algorithm 1. The results for the binomial parameter for N = 10 are plotted in Figure 4.
VI.

CONCLUSION

Optimality of estimators is an always-relevant topic in statistics. Even when optimal estimators are intractable
or impractical, they provide a useful benchmark, and make it possible to show rigorously that some more tractable
estimator is “good enough”. Our main technical contributions in this paper are simple, constructive formulas to
compute Bayes estimators for “fidelity”-type loss functions based on the Bhattacharyya coefficient. This result may
be directly useful for Bayesian machine learning, quantum state estimation, and other inference problems where it
(1) provides a more accurate estimator, and (2) establishes a provable bound on performance.
We find the examples provided in the second half of our paper interesting because they suggest certain qualitative
properties of multinomial estimation. Two of the most widely used estimators for this problem are the MLE and the
posterior mean (which is Bayes-optimal for Bregman divergences). Our results indicate that the Bayes estimator for
Bhattacharya loss interpolates between the MLE and the posterior mean; it hedges away from p = 0 like the mean,
but not as much. Most interestingly, our analysis shows that all the estimators we considered have nearly identical
average risk—and therefore suggests that worrying about finding exact Bayes estimators may be unnecessary.
Superficially, this may appear to undercut our result. Who cares about deriving the Bayes estimator if it isn’t significantly better? But this is exactly the point: our result proves that the posterior mean is “good enough”, at least in this

8
Algorithm 1 Kempthorne’s numerical algorithm for finding the least favorable prior [2].
Input: Number of measurements N > 0.
Input: Support points of initial guess p.
Input: Probability weights of the support points w.
Input: Tolerance tol > 0 (default 10−3 ).
Input: Mixing parameter α (default 0.01).
Output: Least favorable prior { p, w}.
Output: Lower bound on the minimax risk avg risk.
Output: Upper bound on the minimax risk max risk.
function K EMPTHORNE(N, { p, w},tol)
diff ← tol + 1 # make sure we enter the loop
while diff > tol do
{ p, w} ← prior with same number of support points which maximizes the Bayes risk
avg risk ← the maximum value of the Bayes risk for the prior found above
max risk ← global maximum of risk using the Bayes estimator of { p, w}
diff ← |avg risk − max risk|/avg risk
if diff > tol then
Add a new support where the maximum risk is attained
wlength( p) ← α
for each i ≤ length( p) − 1, wi ← wi − α/(length( p) − 1)
end if
end while
return { p, w}, avg risk, max risk
end function

0.5

Average 1 − B 2

0.4

p̂

0.3
0.2

Bayes
Minimax
Mean

0.1
0.0

0

1

2

n

3

4

0.026
0.024
0.022
0.020
0.018
0.016
0.014
0.012
0.010
0.008
0.0

5

Bayes
Minimax
Mean
0.1

0.2

p

0.3

0.4

0.5

FIG. 4: In this figure, we compare three estimators: (1) the minimax estimator, which is the Bayes estimator for the least favorable
prior; (2) the Bayes estimator for the least favorable conjugate prior, given by β ≈ 0.44; and (3) the posterior mean estimator for
the least favorable [with respect to the posterior mean estimator] conjugate prior, given by β ≈ 0.26. On the left, we show the
estimators themselves. On the right, we show their pointwise fidelity (1-risk) as a function of p0 . N = 10 in all cases.

particular case. Furthermore, it shows that although the Bayes and mean estimators are visibly different, their performance is not. This suggests that—at least for Bhattacharya loss functions—a rather wide range of estimators achieve
near-optimal performance. This conclusion is reinforced by the behavior of the minimax (frequentist-optimal) estimators that we construct using our result, where we found that small changes in the precise definition of “minimax”
produced fairly large changes in the “optimal” estimator. And, while we did not fully solve the quantum problem
(by finding Bayes estimators for quantum fidelity), we hope that our partial solution (for commuting states) provides
a stepping stone to a full solution in the future.
Acknowledgments

CF acknowledges funding from the IARPA MQCO program, the ARC via EQuS project number CE11001013, and
by the US Army Research Office grant numbers W911NF-14-1-0098 and W911NF-14-1-0103. Sandia National Laboratories is a multi-mission laboratory managed and operated by Sandia Corporation, a wholly owned subsidiary of

9
Lockheed Martin Corporation, for the U.S. Department of Energy’s National Nuclear Security Administration under
contract DE-AC04-94AL85000.

[1] J. Berger, Statistical Decision Theory and Bayesian Analysis, Springer (1985).
[2] P. J. Kempthorne, Numerical specification of discrete least favorable prior distributions, SIAM Journal on Scientific and Statistical
Computing 8, 171 (1987).
[3] A. Banerjee, Xin Guo and Hui Wang, On the optimality of conditional expectation as a Bregman predictor, IEEE Transactions on
Information Theory 51, 2664 (2005) .
[4] A. Djouadi,O. Snorrason and F. D. Garber, The quality of training sample estimates of the Bhattacharyya coefficient, IEEE Transactions on Pattern Analysis and Machine Intelligence 12, 92 (1990).
[5] W. K. Wootters, Statistical distance and Hilbert space, Physical Review D 23, 357 (1981).
[6] R. Jozsa, Fidelity for mixed quantum states, Journal of Modern Optics 41, 2315 (1994).
[7] C. Ferrie and R. Blume-Kohout, Estimating the bias of a noisy coin, AIP Conference Proceedings 1443, 14 (2012).
[8] E. L. Lehmann and G. Casella, Theory of point estimation, Springer (1998).
[9] J. Aitchison, Goodness of prediction fit, Biometrika 62, 547 (1975).
[10] A. Bhattacharyya, On a measure of divergence between two statistical populations defined by their probability distributions, Bulletin
of the Calcutta Mathematical Society 35, 99 (1943).
[11] S. Boyd, Convex Optimization, Cambridge University Press (2004).
[12] M. A. Nielsen and I. L. Chuang, Quantum computation and quantum information. Cambridge University Press (2010).
[13] L. M. Artiles, R. D. Gill and M. I. Guta, An invitation to quantum tomography, Journal of the Royal Statistical Society: Series B
(Statistical Methodology), 67, 109 (2005).
[14] R. Kueng and C. Ferrie, Near-optimal quantum tomography: estimators and bounds, New Journal of Physics 17, 123013 (2015).
[15] B. S. Clarke and A. R. Barron, Jeffreys’ prior is asymptotically least favorable under entropy risk,Journal of Statistical planning and
Inference 41, 37 (1994).

