3.4 Bayes Estimation

There are two main schools of statistical inference. The frequentist school regards a parameter as an unknown constant. The Bayesian school regards \(\theta \) as a random variable with a prior distribution \(\pi (\theta )\), expressing what is believed about \(\theta \) before the data are seen. The data then update the prior to a posterior distribution, by Bayes’ theorem.

Theorem 3.4.1 (Posterior distribution). If \(\theta \) has prior density \(\pi (\theta )\) and, given \(\theta \), the data \(x = (x_1, \dots , x_n)\) are a random sample from \(f_\theta \), then the conditional density of \(\theta \) given \(x\) is \[\pi (\theta \mid x) = \frac {\pi (\theta )L(\theta )}{\int \pi (t)L(t)\,dt},\qquad L(\theta ) = \prod ^n_{i = 1}f_\theta (x_i) .\]

Proof. The joint density of \((\theta , X)\) is \(\pi (\theta )L(\theta )\), and the marginal density of \(X\) at \(x\) is \(\int \pi (t)L(t)\,dt\). Their ratio is the conditional density. □

In practice one writes \(\pi (\theta \mid x)\propto \pi (\theta )L(\theta )\), that is, posterior is proportional to prior times likelihood, and identifies the normalising constant at the end by recognising the form of the density. Factors of \(L(\theta )\) that do not involve \(\theta \) can be dropped.

Example 3.4.2. A coin with probability \(\theta \) of heads is tossed \(n\) times and gives \(x\) heads. With the prior \(\theta \sim \text {BETA}(10, 10)\), find the posterior distribution of \(\theta \).

Solution. \(L(\theta )\propto \theta ^x(1 - \theta )^{n - x}\) and \(\pi (\theta )\propto \theta ^9(1 - \theta )^9\) on \((0, 1)\), so \[\pi (\theta \mid x)\propto \theta ^{x + 9}(1 - \theta )^{n - x + 9},\qquad 0<\theta <1 .\] This is a beta density, and the constant must be the reciprocal of \(B(x + 10, n - x + 10) = \int _0^1\theta ^{x + 9}(1 - \theta )^{n - x + 9}\,d\theta \), where \(B(\alpha , \beta ) = \frac {\Gamma (\alpha )\Gamma (\beta )}{\Gamma (\alpha + \beta )}\). So \(\theta \mid x\sim \text {BETA}(x + 10, n - x + 10)\).

Definition 3.4.3. A family of prior distributions is conjugate for a model if the posterior distribution is in the same family as the prior for every data set.

Example 3.4.2 shows that the beta family is conjugate for the binomial. The prior parameters act like extra data: a \(\text {BETA}(10, 10)\) prior behaves like \(20\) earlier tosses with \(10\) heads.

Example 3.4.4. Let \(X_1, \dots , X_n\) be a random sample from \(f_\theta (x) = \theta x^{\theta - 1}\), \(0<x<1\), \(\theta >0\). Show that the prior \(\theta \sim \text {GAM}\left (a + 1, \frac 1b\right )\), with density proportional to \(\theta ^ae^{-b\theta }\), is conjugate.

Solution. \(L(\theta ) = \theta ^n\left (\prod x_i\right )^{\theta - 1}\propto \theta ^n e^{\theta \sum \ln x_i}\), so \[\pi (\theta \mid x)\propto \theta ^ae^{-b\theta }\,\theta ^ne^{\theta \sum \ln x_i} = \theta ^{a + n}\exp \left \{-\theta \left (b - \sum \ln x_i\right )\right \}.\] Since every \(\ln x_i<0\), the coefficient \(b - \sum \ln x_i\) is positive, and the posterior is \[\theta \mid x\sim \text {GAM}\left (a + n + 1, \frac {1}{b - \sum \ln x_i}\right ),\] again a gamma distribution.

Bayes point estimation

A point estimate can be obtained from the posterior distribution by choosing a loss function and minimising the expected loss under the posterior. Which estimator results depends on the loss.

Theorem 3.4.5. For squared error loss \(L(\theta , t) = (t - \theta )^2\), the Bayes estimator, the value \(t\) that minimises the posterior expected loss \(E\left [(t - \theta )^2\mid x\right ]\), is the posterior mean \[\tilde \theta (x) = E(\theta \mid x) = \int \theta \,\pi (\theta \mid x)\,d\theta .\]

Proof. Write \(m = E(\theta \mid x)\). Then \[E\left [(t - \theta )^2\mid x\right ] = E\left [(t - m + m - \theta )^2\mid x\right ] = (t - m)^2 + \var (\theta \mid x),\] since the cross term \(2(t - m)E(m - \theta \mid x)\) is zero. This is smallest at \(t = m\). □

Theorem 3.4.6. For absolute error loss \(L(\theta , t) = |t - \theta |\), the Bayes estimator is any median of the posterior distribution.

Proof. Let \(F\) be the posterior distribution function and \(g(t) = E\left [|t - \theta |\mid x\right ] = \int _{-\infty }^t(t - \theta )\,dF(\theta ) + \int _t^\infty (\theta - t)\,dF(\theta )\). For a continuous posterior, differentiating gives \(g'(t) = F(t) - (1 - F(t)) = 2F(t) - 1\), which is negative for \(F(t)<\frac 12\) and positive for \(F(t)>\frac 12\). So \(g\) is minimised where \(F(t) = \frac 12\). □

Example 3.4.7. In example 3.4.4, find the Bayes estimator of \(\theta \) under squared error loss, and compare it with the ML estimator.

Solution. A \(\text {GAM}(\alpha , \beta )\) distribution has mean \(\alpha \beta \), so the posterior mean is \[\tilde \theta = \frac {a + n + 1}{b - \sum \ln X_i} .\] With \(T = -\sum \ln X_i\), the ML estimator is \(\frac nT\) (example 2.5.5) and the Bayes estimator is \(\frac {n + a + 1}{T + b}\). The prior adds \(a + 1\) to the numerator and \(b\) to the denominator, and its influence fades as \(n\) and \(T\) grow.

Example 3.4.8. In example 3.4.2, suppose \(n = 20\) tosses give \(x = 14\) heads. Find the Bayes estimates of \(\theta \) under squared and absolute error loss, and compare them with the ML estimate.

Solution. The posterior is \(\text {BETA}(24, 16)\). Its mean is \(\frac {24}{24 + 16} = 0.600\), and its median, found numerically, is \(0.602\) (for a beta distribution with both parameters above \(1\) the median is close to \(\frac {\alpha - 1/3}{\alpha + \beta - 2/3} = 0.602\)). The ML estimate is \(\frac {14}{20} = 0.70\). The posterior mean can be written as \[\frac {24}{40} = \frac {20}{40}\cdot \frac {14}{20} + \frac {20}{40}\cdot \frac {10}{20},\] a weighted average of the sample proportion \(0.7\) and the prior mean \(0.5\), with weights proportional to the sample size \(20\) and the prior’s equivalent sample size \(20\). This shrinkage towards the prior mean is typical of Bayes estimators.

Problem 3.4.1. Let \(X_1, \, X_2, \, \cdots \, , \, X_n\) be a random sample from the \(POI(\theta )\) distribution.

(a): Show that the prior distribution \[\theta \, \thicksim \, GAM\left (b + 1, \frac {1}{a}\right )\] is a conjugate prior.
(b): using the conjugate prior in (a) find the Bayes estimator of \(\theta \) for a squared error loss function.

Show solution

Solution.

(a) Conjugacy

Write \(T=\sum _{i=1}^{n}X_i\). The likelihood and the GAM\(\left (b+1,1/a\right )\) prior are \[L(\theta ) \propto e^{-n\theta }\theta ^{T}, \qquad \pi (\theta ) \propto \theta ^{b}e^{-a\theta },\quad \theta >0 .\] Multiplying, \[\pi \left (\theta \mid x\right ) \propto \theta ^{\,T+b}\,e^{-(n+a)\theta },\] which is the kernel of a GAM\(\left (T+b+1,\ \frac {1}{n+a}\right )\) density. The posterior belongs to the same gamma family as the prior, so the family is conjugate for the Poisson likelihood.

(b) The Bayes estimator

Under squared error loss the Bayes estimator is the posterior mean. A GAM\((\kappa ,\lambda )\) variable has mean \(\kappa \lambda \), so \[\widetilde {\theta } = \frac {T+b+1}{n+a} = \frac {\sum _{i=1}^{n}X_i + b + 1}{n+a} .\]

Note. Write it as a weighted average: \[\widetilde {\theta } = \frac {n}{n+a}\,\overline {X} + \frac {a}{n+a}\cdot \frac {b+1}{a},\] a compromise between the sample mean and the prior mean \((b+1)/a\), with weights \(n\) and \(a\). So \(a\) acts as a prior sample size: the prior carries as much influence as \(a\) observations, and its influence vanishes as \(n\) grows.

Problem 3.4.2. Let \(X_1, \, X_2, \, \cdots \, , \, X_n\) be a random sample from the \(X\thicksim N(0, \sigma ^2)\) distribution where \(\sigma ^2\) is known. Assuming the prior distribution is \(\theta \, \thicksim \, N\left (\frac {b}{a}, \frac {\sigma ^2}{a}\right ),\) fined the Bayes estimator using a squared error loss function.

Show solution

Solution. The model must be \(X_i \sim N\left (\theta ,\sigma ^{2}\right )\) with \(\sigma ^{2}\) known — as printed the mean is zero and \(\theta \) does not appear in the likelihood at all, so there would be nothing to estimate.

With \(T=\sum _i X_i\) and prior \(\theta \sim N\left (b/a,\ \sigma ^{2}/a\right )\), \[\pi \left (\theta \mid x\right ) \propto \exp \left \{-\frac {1}{2\sigma ^{2}}\sum _i\left (x_i-\theta \right )^{2}\right \} \exp \left \{-\frac {a}{2\sigma ^{2}}\left (\theta -\frac {b}{a}\right )^{2}\right \}.\] Both exponents are quadratic in \(\theta \), so the posterior is normal and it is enough to collect the coefficients. Discarding terms free of \(\theta \), the exponent is \[-\frac {1}{2\sigma ^{2}}\left [(n+a)\theta ^{2} - 2\theta \left (T+b\right )\right ] = -\frac {n+a}{2\sigma ^{2}}\left [\theta - \frac {T+b}{n+a}\right ]^{2} + \text {const} .\] Hence \[\theta \mid x \ \sim \ N\!\left (\frac {T+b}{n+a},\ \frac {\sigma ^{2}}{n+a}\right ),\] and under squared error loss the Bayes estimator is the posterior mean \[\widetilde {\theta } = \frac {\sum _{i=1}^{n}X_i + b}{n+a} .\]

As in the Poisson problem this is a weighted average of \(\overline {X}\) and the prior mean \(b/a\), with weights \(n\) and \(a\); the prior is worth \(a\) observations, and the posterior variance \(\sigma ^{2}/(n+a)\) says the same thing.

Problem 3.4.3. Let \(X_1, \, X_2, \, \cdots \, , \, X_n\) be a random sample from the \(POI(\theta )\) distribution. Assuming the prior distribution of \(\theta \) is \(\theta \thicksim BETA(1/2,1)\), find the Bayes estimator of \(\theta \) using the squared error loss function.

Show solution

Solution. The BETA\((1/2,1)\) prior is \(\pi (\theta )\propto \theta ^{-1/2}\) on \((0,1)\), so with \(T=\sum _i X_i\), \[\pi \left (\theta \mid x\right ) \propto \theta ^{\,T-1/2}e^{-n\theta }, \qquad 0<\theta <1 .\] This is a gamma kernel truncated to the unit interval, not a gamma density, so the normalising constant is an incomplete gamma function. Using \[\int _0^{1}\theta ^{\,s-1}e^{-n\theta }\,d\theta = n^{-s}\,\gamma (s,n),\] where \(\gamma \) is the lower incomplete gamma function, the posterior mean is \[\widetilde {\theta } = \frac {\int _0^{1}\theta ^{\,T+1/2}e^{-n\theta }d\theta }{\int _0^{1}\theta ^{\,T-1/2}e^{-n\theta }d\theta } = \frac {n^{-(T+3/2)}\gamma \left (T+\tfrac 32,\,n\right )}{n^{-(T+1/2)}\gamma \left (T+\tfrac 12,\,n\right )} = \frac {\gamma \left (T+\tfrac 32,\,n\right )}{n\,\gamma \left (T+\tfrac 12,\,n\right )} .\] This is the Bayes estimator under squared error loss. It has no expression in elementary functions and is evaluated numerically.

Remark. The prior deserves comment. A BETA\((1/2,1)\) distribution lives on \((0,1)\), so it asserts with certainty that the Poisson rate is below one. That is a strong claim, and it is what destroys conjugacy — the gamma prior of the previous problem is supported on all of \((0,\infty )\) and returns a gamma posterior. If the restriction is genuinely believed the answer above is correct; if it was chosen for convenience it is the wrong prior, and it is worth noticing that the inconvenience of the answer is the model telling us so.

Questions on this section

Stuck on something here? Ask below and it stays attached to this topic.