2.3.1 The Likelihood Principle

Let us consider an i.i.d. sample x = (x₁, …, x_n) from an infinite population, characterized by a distribution density function f(x, θ), where x is the vector of data and θ is a (maybe, vector) parameter. Recall from basic statistics course that the likelihood function is defined as

(2.13)

In the context of incomplete data, which is typical for many risk management and reliability applications, including some of those which will be used in the following chapters, we may observe exact values x = (x₁, …, x_m), m < n and just know the lower limits for the other cases: x_j > c_j, j = m + 1, …, n. This corresponds to the so-called right censoring of the dataset. In this case we can modify the likelihood formula:

(2.14)

where S(c, θ) = P(X > c ; θ) is the survival function corresponding to the distribution with density f(x, θ).

The likelihood principle is a practically undisputed basic principle of mathematical statistics. According to the likelihood principle, all information regarding the unknown parameter, which could be found in data, is contained in the likelihood function.

2.3.2 Review of Classical Procedures

Let us review the main definitions and procedures of hypothesis testing for the simplest case of two parametric hypotheses. We believe that there exist two exhaustive and mutually exclusive possibilities regarding the parameter θ of the distribution f(x, θ). One is designated as the null hypothesis and is denoted by H₀: θ = θ₀. The other is described as the alternative H₁: θ = θ₁. A test statistic (a function of the sample) T(x) and critical region R_T⊂Range(T) are used to define a test, which is a basic binary decision rule:

“Reject H₀” if T(x) ∈ R_T or “Do not reject H₀” if T(x)∉R_T.

Notice that although in the formal context of two simple hypotheses “Do not reject H₀” seems to be equivalent to “Accept H₀” or “Reject H₁,” in common applications we find the first of the three statements the most appropriate. Null hypothesis is also known as status-quo or zero-effect hypothesis requiring no immediate action if not rejected. Typical examples of a null hypothesis would be: “new drug is no good,” “surgery is not required,” “new technology is exactly as effective as the old one.” Rejecting null is something which requires an action: “approve a new drug,” “proceed with the surgery,” “implement the new technology” and is associated with additional cost or risk. The burden of proof is higher for rejecting the null. Not rejecting the null might call for more testing to be performed and more proofs to be provided before we take the action. It does not necessarily mean accepting the null or rejecting the alternative. In the context of classical hypotheses testing, one usually speaks of statistical significance of the results, if the null hypothesis can be rejected.

We can judge the quality of the tests using definitions of two errors which we want to be as small as possible:

• Type I error (reject true null) with probability α = P(T(x) ∈ R_T∣H₀)
• Type II error (do not reject false null) with probability β = P(T(x)∉R_T∣H₁)

Quantity 1 − β is also known as the power of the test.

In general, type I error is the one which is more important in the context of the study. Therefore, a common approach is to choose the largest admissible probability of type I error, which is called significance level, and then find tests which will make type II error as small as possible. Neyman–Pearson theory of hypothesis testing suggests that in a very general setting the most powerful tests (those minimizing type II error for type I error not exceeding the significance level) should be based on the likelihood ratio:

(2.15)

If we use the likelihood ratio as the test statistic, the only thing we need to specify the test procedure is the critical region, which will take the form R_Λ = {x: Λ(x) > c_α}. Here c_α is the critical value of the test, which can be determined by α, if the distribution of the likelihood ratio under the condition of null hypothesis is known. In practical applications, we usually look for a simpler test statistic T(x) equivalent to the likelihood ratio in the sense that R_T = {x: T(x) > k_α} = R_Λ = {x: Λ(x) > c_α}.

Ten Coin Tosses

A random variable X has Bernoulli distribution (see Section 1.2) if it may assume two values only: 0 or 1, and is described by a discrete probability function f(x, θ) = θ^x(1 − θ)^{1 − x}, where the value of the parameter θ = P(X = 1) is unknown. Consider a random i.i.d. sample x = (x₁, …, x_n) from this distribution. Write down the likelihood function

(2.16)

Let us perform a series of ten trials: tosses of the same coin. Suppose we record “heads” seven times. This series is treated as an i.i.d. sample x = (x₁, …, x_n), n = 10, where the underlying variable X has Bernoulli distribution and θ is unknown probability of “heads” (success). We will consider null hypothesis H₀: θ = θ₀ = 0.5 (fair coin) versus a simple alternative H₁: θ = θ₁ = 0.7 (coin biased in a special way). How can we construct the likelihood ratio test?

In our case,

(2.17)

so

(2.18)

The exact formula with logarithms on the right-hand side of the last expression is irrelevant as soon as we can see that for some k_α, and the equivalent form of the likelihood ratio test is:

“Reject H₀” if or “Do not reject H₀” if .

Picking up an “industry standard” α = 0.05, from the table of binomial distribution we can see that, P(∑X_i > 7∣H₀) = 0.055, P(∑X_i > 8∣H₀) = 0.011, and the rejection region consists of so in our example with seven heads we do not reject the null hypothesis that the coin is fair. Notice also that our final verdict does not depend on the alternative if instead of θ₁ = 0.7 we choose another value larger than 0.5. The test might have less power, but its significance level stays the same.

2.3.3 Bayesian Hypotheses Testing

Consider the same general setting of two simple hypotheses. Incorporating prior information in this case means assigning prior probabilities to the null hypothesis P(H₀) and the alternative P(H₁) = 1 − P(H₀).

Let us also consider the observed data x and the likelihoods L(x, θ₀) and L(x, θ₁). Since the two hypotheses are mutually exclusive and exhaustive, the Law of Total Probability suggests to “integrate out” the parameter and obtain a parameter-free weighted likelihood, where the prior probabilities of the hypotheses play the role of weights:

As a function of data, m(x) is usually called the marginal distribution. For given data it is a constant, and as we will see very soon, it plays very little role in the following analysis. However, we will need it during the derivation of the Bayesian test. According to the definition of conditional probability or a simplified version of the Bayes theorem,

and the ratio of these two expressions, after m(x) cancels out, is

(2.19)

The ratio on the left-hand side is known as posterior odds. The first ratio on the right-hand side is the likelihood ratio, which we have encountered while developing classical testing procedures. The second is the ratio of prior probabilities or prior odds. When the prior odds equal one, the posterior odds has the Bayesian name of the Bayes factor, see [15]. In this case of two simple hypotheses, the Bayes factor coincides with the likelihood ratio. In general, both the Bayes factor and the prior odds will be needed to evaluate the posterior odds.

Finally, in order to decide whether to reject the null hypothesis or not, without any prior information (that means, setting prior odds to one, or in other words, assuming P(H₀) = P(H₁) = 0.5), we have to establish a reasonable benchmark for the posterior odds. We can consider the posterior odds 3:1 to provide “substantial” evidence against the null and 10:1 to provide “strong” evidence. Jeffreys [14] used these benchmarks for the Bayes factors, see also [15]. So for substantial evidence against the null, the Bayesian test takes the form:

“Reject H₀” if or “Do not reject H₀” if

For strong evidence,

“Reject H₀” if or “Do not reject H₀” if

Ten Coin Tosses

Returning to the ten coin toss experiment in Bayesian setting, we need to specify the priors. Let us say that our prior belief in fairness of the coin can be expressed by P(H₀) = 0.5 or in the prior odds of 1:1. From (2.17) we may calculate the Bayes factor as

bringing about the posterior odds of 2.28:1 which gives less than substantial (and much less than strong) evidence against the null hypothesis. Notice that the classical procedure introduced in Subsection 2.3.2 did not require the computation of the likelihood ratio, but it called for the computation of binomial probabilities.

To extend the analogy with the classical procedure, assume that for the same null and alternative and choice of priors we observe a different result: instead of seven, we have nine “heads” in a series of ten coin tosses. Then the Bayes factor is

which provides a strong evidence against the null. Noticing that lies in the critical region for α = 0.05, we may conclude that in this case both classical and Bayesian approaches to hypotheses testing lead us to similar results.

2.4.1 Review of Classical Procedures

Let us recall some concepts of classical theory of statistical estimation. Consider a random sample x = (x₁, …, x_n) from a certain population with p.d.f. f(x, θ). A statistic is defined as any function of the random sample. A statistic used as an approximate value of unknown parameter of the distribution (population) can be considered an estimator of the parameter. All estimators are random variables. We have to distinguish an estimate as a specific numeric value used to approximate the unknown parameter from an estimator which is a rule leading to different numerical estimates for different samples from the same distribution (population). We will try not to mix up these two terms, though it might be confusing at times.

From the classical standpoint, desirable properties of estimators are: consistency, unbiasedness, and efficiency. An estimator is consistent if it converges in a certain sense to the true parameter value when the sample size grows to infinity. An estimator is unbiased if its expected value (mean) coincides with the true value of the estimated parameter. An estimator is called the most efficient within a certain class if for the given class its mean square error is the smallest.

2.4.2 Maximum Likelihood Estimation

There exist several popular methods to obtain parametric estimators in a general setting: the method of moments, the method of least squares, etc. However the method of maximum likelihood is probably the most important for classical statistics. A maximum likelihood estimator (MLE) for unknown k-dimensional vector parameter θ is a statistic, which maximizes the likelihood function defined in (2.13). When the sample sizes are large, or ideally, are allowed to grow to infinity, MLEs have all nice properties: they are consistent, asymptotically unbiased, asymptotically efficient, and even asymptotically normal (their distribution is close to normal). However for small sample sizes MLEs do not always behave that well, and may be substantially inferior to the estimators obtained by other methods.

Let us recall the algorithm used to obtain MLE in regular cases.

• Write down the likelihood function (2.13).
• Calculate logarithmic likelihood or log-likelihood l(x₁, …, x_n, θ) = ln L (a technical step simplifying further derivation).
• Calculate partial derivatives of the log-likelihood with respect to parameters.
• Find the critical points of the logarithmic likelihood, that is, the solutions of the system of equations
  (2.20)
• Check if the critical point found is indeed the absolute maximum point (not a local maximum, nor a minimum or a saddle point). The value of absolute maximum obtained herewith is the MLE.

This algorithm cannot be used in nonregular cases (e.g., if the partial derivatives are not defined). Sometimes it does not lead to a unique solution. Notice also that the exact solution of system (2.20) is available only in some simple textbook cases. In more complex situations one has to use numerical methods of optimization to find an approximate solution. Let us recall two simple examples illustrating both analytical and numerical methods.

Ten Coin Tosses

Returning to the ten coin toss experiment, we will once again use the likelihood function defined in (2.16). Then we will calculate its logarithm, the log-likelihood

(2.21)

It is easy to find the derivative w.r.t θ and solve the following equation

(2.22)

The solution is not unexpected: . To check if this critical point is a maximum, find the second derivative

(2.23)

Since for Bernoulli distributed samples , this second derivative is never positive. Therefore, the critical point we have found is a local maximum point. As soon as it is the only critical point, it is also the absolute maximum. Thus is the MLE for the parameter of Bernoulli distribution.

We perform a series of 10 trials (coin tosses): n = 10. Suppose we record “heads” seven times (x = 7). So .

Life Length

Many random variables in insurance or reliability applications representing “life length” or “time to failure” can be conveniently described by gamma distribution with density function

(2.24)

where the values of two parameters α (shape) and λ (rate or reciprocal scale) are unknown, see also Section 1.4. Gamma class of distribution was introduced in Chapter 1. The variety of forms of gamma densities for different combinations of parameter values (see Figure 1.17) makes this class very attractive in modeling survival.

Suppose we obtain an i.i.d. sample from gamma distribution (x₁, …, x_n). Let us write down the log-likelihood function

(2.25)

As we can see, the approach of the previous example will not work that easily: differentiation in both α and λ is way too difficult to allow for a simple analytical solution. However, numerical optimization procedures might help (see Exercises at the end of the chapter).

An alternative to MLE is known as the method of moments, which is a particular case of a more general plug-in principle. It uses the fact that if there exists a convenient formula expressing the unknown parameters of a distribution through its moments (say, the mean and the variance), then plugging sample moments in instead of theoretical moments often brings about consistent (not necessarily unbiased or efficient) estimators.

From (1.19) we know that E(X) = α/λ and Var(X) = α/λ², therefore

after plug-in leading to the method of moment estimators:

Method of moments estimators sometimes are attractively simple, but are lacking certain large sample optimality properties of MLE. Returning to 10 coins example, we can easily see that due to E(X) = θ, the method of moments estimator of θ coincides with the MLE.

2.4.3 Bayesian Approach to Parametric Estimation

Let us assume that all relevant prior information regarding parameter θ is summarized in π(θ), which is the p.d.f. of the prior distribution. In order to obtain the posterior distribution, we can use continuous version of the Bayes Theorem, the most important formula in this chapter!

(2.26)

The left-hand term of (2.26) is the p.d.f. of the posterior distribution. The right-hand term is just a product of prior and likelihood. Symbol ∝ denotes proportionality of its left and right hand sides or equivalence of those to within a constant multiple. The constant multiple in question is the integral in the denominator of the central term of (2.26). This integral depends on x only and does not depend on the parameter which is integrated out. It represents the marginal distribution of the data, and is an absolute constant when data are provided. The bad news is: In most applications this constant is hard to determine. The good news is: In most cases it is not necessary. There exist some ways around.

For now let us suppose that this integral is not a problem, the posterior distribution is defined completely, and it has a nice analytical form. This is possible in a limited number of simple examples. The posterior distribution is what we need to provide Bayesian inference. What information can we get out of the posterior distribution?

Point Estimation

The mean, median, or mode of the posterior distribution (in the sequel: the posterior mean, the posterior median, and the posterior mode) can play the role of point estimators of the parameter:

(2.27)

(2.28)

(2.29)

The first one (by far, the most popular choice) minimizes posterior Bayes risk which is defined as the mean square deviation of the estimator from the parameter averaged over all parameter values, while the second minimizes the absolute deviation of the estimator. Sometimes, the posterior mode is the easiest to use.

Interval Estimation

Having chosen the level 1 − α in advance, we are looking for an interval (θ_l, θ_u) such that

(2.30)

Such an interval is called a credible interval. Unlike classical confidence intervals, Bayesian credible intervals allow for a simple probabilistic interpretation as intervals of given posterior probability. We can always choose a symmetric interval, though it is not the only possible choice:

(2.31)

2.5.1 Classical (Frequentist) Approach

Classical approach, in our view, found its most logical representation in the ideas and works of a prominent British statistician Sir Ronald Aylmer Fisher (1890–1962) (Figure 2.4). As the portrait reveals, Fisher was an avid smoker and did some research on the effect of smoking tobacco on human health. What is more important, he made a great contribution to the development of Statistics in the twentieth century. He worked in multiple applications and obtained many deep theoretical results. He also contributed to the formulation of basic principles of what we call here the classical approach. He also happened to be very anti-Bayesian.

Let us recall the main principles of classical statistical analysis:

• Parameters (numerical characteristics of populations and processes) are fixed unknown constants.
• Probabilities are always interpreted as the limits of relative frequencies in unlimited series of repetitive experiments.
• Statistical methods are judged by their large sample behavior.

A typical method of statistical estimation is the maximum likelihood method. Hypotheses testing is based on the likelihood ratio. The likelihood function is treated as a function of data and is not considered to describe the distribution of the parameter.

However, in order to cast some doubt in practical infallibility of this approach, we may use a classical experiment performed by Fisher and discussed in his book [10] to illustrate Fisher’s exact test. In our description of this experiment we will follow [20], though this story being told by many statisticians has many different versions.

2.5.2 Lady Tasting Tea

At a summer tea party in Cambridge, England, a lady (some sources identify her as Dr. Muriel Bristol) states that tea poured into milk tastes differently than that of milk poured into tea. The lady claims to be able to tell whether the tea or the milk was added first to a cup. Her notion is taken sceptically by the scientific minds of the group. But one guest, happening to be Ronald Aylmer Fisher, proposes to scientifically test the lady’s hypothesis. Fisher proposes to give her eight cups, four of each variety, in random order. One could then ask what the probability was for her getting the number she got correct, but just by chance (the null hypothesis). The experiment provides the lady with eight randomly ordered cups of tea, four prepared by first adding milk, four prepared by first adding the tea. She is to select the four cups prepared by one method. She was reported [22] to prove correct all eight times. Using permutation formula, one can suggest that under the null hypothesis her chances to guess correctly, or randomly select one out of 8!/(4!(8 − 4)!) = 70 possible permutations of eight cups are about 0.014, which is the p-value of Fisher’s exact test.

Putting aside interesting connotations of this problem for experimental design and randomization (the latter to be discussed in the next chapter), we can consider a slightly different setup, where each of the eight cups is randomly and independently prepared in one of the two possible fashions (the milk first or the tea first) according to the results of a fair coin toss. Then the chances of the lady to guess correctly all eight times become 1/2⁸ = 1/256, which is approximately 0.004. If we increase the number of cups from 8 to 10, then in the case of all 10 correct guesses the p-value or the significance level of the test goes down to 0.001. However, what exactly does this p-value mean?

Let us cite a letter from an outstanding American Bayesian statistician L. J. Savage (Figure 2.5) to a well-known colleague, biostatistician J. Cornfield, where this experiment (discussing a version with 10 cups) is used to cast some doubt on the practicality of using classical procedures based solely on likelihood and suggesting some role of the context of the experiment and thus the prior information.

“There is a sense in which a given likelihood ratio does have the same sense from one problem to another, at least from the Bayesian point of view. It is always the factor by which the posterior odds differ from the prior odds. However, all schools of thought are properly agreed that the critical likelihood ratio will vary from one application to another. The Neyman-Pearson school expresses this by saying that the choice of the likelihood ratio is subjective and should be made by the user of the data for his particular purpose. As a Bayesian, I would say that a critical likelihood ratio for a dichotomous decision about a simple dichotomy depends (in an evident way) on the loss associated with the decision and on the prior probabilities associated with the dichotomy.

The fact that the reaction to an experiment depends on the content of the experiment as well as on its mathematical structure and whatever economic issues might be involved seems to me to be brought out by the following triplet of examples that occurred to me the first time I taught statistics, when I still had a completely orthodox orientation. Imagine the following three experiments: 1. Fisher’s lady has correctly dealt with ten pairs of cups of tea. Many readers will recognize this example to be from R. A. Fisher’s famous lady tasting tea experiment. 2. The professor of 18-th century musicology at the University of Vienna has correctly decided for each of 10 pairs of pages of music which was written by Mozart and which by Haydn. 3. A drunk in a parlor has succeeded 10 times in correctly calling a coin secretly tossed by you. These three experiments all have the same mathematical structure and the same high significance level. Can there, however, be any question that your reaction to them is justifiably different? My own would be: 1. I am still skeptical of the lady’s claim, but her success in her experiment has definitely opened my mind. 2. I would originally have expected the musicologist to make this discrimination; I would even expect some success in making it myself; he, an expert in the matter, felt sure that he could make it. His success in 10 correct trials confirms my original judgment and leaves no practical doubt that he would be correct in substantially more than half of future trials, though I would not be surprised if he made occasional errors. 3. My original belief in clairvoyance was academic, if not utterly nonexistent. I do not even believe that the trial was conducted in such a way that trickery is a plausible hypothesis, and feel sure that the drunk simply had an unusual run of luck. Of course, these tests are not simple dichotomies, but I think you will find them germane to your question” [21].

2.5.3 Bayes Theorem

History of the Bayesian approach in statistics traditionally goes back to the times of Reverend Thomas Bayes (Figure 2.6). One can argue though that Laplace and other prominent contemporaries did a lot of what we now call Bayesian reasoning. Not much is known for sure about Bayes including his date of birth. Even his portrait included in this book, according to modern studies, is no longer believed to be his.

Thomas Bayes in his posthumous treatise formulated the statement known now as the Bayes Theorem. This theorem indicates how one should calculate inverse conditional probabilities of hypotheses after the experiment related to these hypotheses. This theorem can be rigorously formulated for the discrete case

(2.32)

Here B_i are hypotheses, and A is some event, which is influenced by the hypotheses, thus the post-event probabilities of the hypotheses are also influenced by the event. This formula in a simple form was used in the Bayesian treatment of two simple hypotheses in Section 2.3. In estimation problems we are more likely to use the continuous version of the theorem which was introduced in (2.26).

Bayes’ work and ideas were never completely forgotten, but they also were out of the mainstream for a very long time. The revival of Bayesian approach is probably due to two factors. The first factor is the recent expansion of new applications of statistics (insurance, biology, medicine) characterized by relatively small sample sizes (limited data) and ample prior information. The second is the development of computational techniques which substantially extends the practical applicability of Bayesian methods.

It would probably be most accurate to name another outstanding British scientist Sir Harold Jeffreys (1891–1989) (Figure 2.7) as the founder of the modern Bayesian statistics. One of his main contributions was his seminal book [14] containing the foundations of Bayesian statistics. However many other great statisticians, including Bruno De Finetti, the famous proponent of subjective probability; David Lindley; already cited Jimmie Savage; and Arnold Zellner (Figure 2.8) should be honorably mentioned among those who brought attention back to this eighteenth century approach.

2.5.4 Main Principles of the Bayesian Approach

Let us summarize the main principles of Bayesian statistical analysis keeping close to the treatment of Bolstad [7].

• True values of parameters are unknown, are not likely to become known, and therefore will be considered random.
• Prior information regarding the parameters before the experiment is allowed and is mathematically expressed in the form of the prior distribution of the random parameter.
• All probabilistic statements expressed by the researcher regarding the parameters express the “degree of rational belief” and are subjective.
• Likelihood function is the only source of information regarding the parameter which is contained in data (weak likelihood principle).
• We do not consider the feasibility of obtaining new data and do not consider the data which could theoretically be obtained.
• Probability rules and laws are applied directly to the statistical inference.
• Our information regarding the parameters is updated using the data.
• Posterior distribution incorporates the prior information and the information from data in the only mathematically consistent way: using the Bayes formula.

The main advantages of Bayesian approach are well known and are discussed in different terms and at different levels in such books as [1, 4, 7, 11], and [16]. From decision theoretic standpoint, we can prove that all admissible estimation procedures can be obtained using Bayesian approach. Sometimes, Bayesian procedures outperform classical procedures even using classical criteria of performance. Bayesian estimation can be substantially more reliable for small samples, while it will lead to the same results as MLE when the sample size grows to infinity. Bayesian procedures behave better than classical ones when there exists substantial amount of prior information which should not be neglected. However, even with very weak prior information, Bayesian approach sometimes performs better than the classical approach. One explanation to this phenomenon is a possible difference in numerical stability between mathematical procedures of optimization (MLE) and integration (Bayesian estimators).

The development of Bayesian procedures in the past faced a very important obstacle. The problem may sound technical, but it was grave. The continuous version of the Bayes theorem (2.26) allows one to determine the posterior density to within a multiplicative constant. Combining (2.26) and (2.27) we can write down the following formula for the most popular Bayes estimator: the posterior mean.

(2.33)

Unfortunately, exact analytical solution for these integrals is rarely available, although we will consider some nice exclusions in Section 2.7. Numerical integration is possible but is far from trivial especially for higher dimensions of parametric space (multiple parameters). That is why Bayesian statistics was stuck in a rut for a long time, being restricted to handling very special cases and not being able to spread the wings.

A “Bayesian revolution” took place in the middle of the twentieth century. What Bayes could not dream of, and what Jeffreys could not carry out, was made possible by the development of the toolkit widely known as MCMC—Markov chain Monte Carlo methods. International Society for Bayesian Analysis (ISBA) was founded in 1992 as the result of increased interest to Bayesian approach, and since then Bayes methods have been growing in popularity even more. Modern Bayesian statistics provides ample tools for such integration developed since the 1950s and still under active development. Chapter 4 of the book is dedicated to some of these methods.

However, the first question we have to ask before we face the problems with evaluating the posterior in (2.33) is: How do we specify the prior?

2.6.1 Subjective Priors

Consider a particular problem with a substantial amount of available prior information: Before the experiment we already know something important about the parameter of interest. This knowledge can be expressed in terms of the parameter’s probability distribution. If probabilities are understood as the degrees of rational belief, no wonder that different people may have different subjective beliefs which will lead to different subjective distributions. Two medical doctors read a patient’s chart and make different observations. Two investors interpret long-term market information in two different ways. They will have different priors.

Assume that both doctors perform the same necessary tests and both investors observe the same results of the recent trades, meaning both parties get the same data and agree on the likelihood. Their final conclusions may still differ because of their priors being different.

This difference in the final decisions due to the difference of priors is both the curse and the blessing of Bayesian statistics. If the use of prior information is allowed, subjectivity is inevitably incorporated in the final decision. Therefore the best thing to do is to develop the rules which will translate subjective prior beliefs into a properly specified prior distribution.

Imagine you are a statistician hired to help with an applied problem: You have to determine how likely a given coin is to fall heads up. One should understand that when we are talking “coin tosses” and “heads,” we do not necessarily have actual coins in mind. We might model binary random events with probability of success being represented by a “head,” and probability of failure represented by a “tail” in a clinical trial setting, as well as while recommending medical treatment or investment decisions.

Statisticians rarely work on applied problems on their own; they become team members along with experts in the object field. An expert on coins joins you on the team and you two have to state the problem statistically and develop an experiment together. Both of you have agreed that available time and resources allow you to perform 10 independent identical tosses of the coin. The parameter of interest is probability of heads θ and the experimental data X = (X₁, …, X_n) are results of tosses, n = 10, P(X_i = 1) = θ, P(X_i = 0) = 1 − θ, X_i being the random number of heads on each toss. This problem was considered above as an illustration of maximum likelihood estimation.

However, now you are allowed to observe the coin prior to the experiment and summarize the prior information. The coin expert can provide the knowledge and you will assist with putting this knowledge into the form of prior distribution on θ. The expert agrees that chances of the coin landing on its edge or disappearing in thin air during any toss are negligible, so that any outcomes other than heads or tails should not be taken into account. He also observes the coin and it looks symmetric. It is likely that the coin is fair, but the past experience tells the expert that some coins are biased either toward the heads or toward the tails. You ask him to quantify these considerations in terms of average expectations. Suppose that the expert concludes: on the average, θ is about 0.5, but this guess is expected to be off on the average by 0.1. You interpret this expert opinion in terms of two numbers: mean value E(θ) = 0.5 and variance Var(θ) = [st.dev.(θ)]² = 0.01.

Now it is your time to specify the prior. The choice of beta family of distributions introduced in Section 1.5 for random θ is logical. It is a wide class including both symmetric and skewed, flat and peaked distributions: θ ∼ Beta(α, β), if its p.d.f. according to (1.25) is

for 0 < θ < 1, α ≥ 0, β ≥ 0. The cluster of gamma functions on the right-hand side does not depend on θ thus it is irrelevant for further analysis so we can write

(2.34)

Several Beta p.d.f.s for various parameter values are depicted in Figure 2.9, see also Figure 1.11.

One of the representatives of the beta family will fit the requirements of E(θ) = 0.5 and Var(θ) = 0.01. Using formulas from Chapter 1 (1.26)

we can solve for α and β:

and

Then the expert’s opinion can be paraphrased as θ ∼ Beta(12, 12). This is our subjective prior. We will use this example in future because it works so well numerically.

Ten Coin Tosses

To complete our analysis and obtain the Bayes estimate for θ, consider an experiment with 7 heads in 10 tosses, . Using (2.16), we obtain the likelihood

(2.35)

and using (2.26), the posterior

The posterior is also a beta distribution shown along with the prior in Figure 2.10, so using the first part of formula (2.35) with new values α = 19 and β = 15, we obtain the numerical value of the Bayes estimate for θ

(2.36)

while the second part of (1.26) provides an error estimate

(2.37)

2.6.2 Objective Priors

The main benefit of the Bayesian approach to statistics is the ability to incorporate prior information. Suppose no prior information is available, and even object field experts cannot help. What could be the rationale for applying Bayesian approach instead of classical in such cases anyway? Why should we do it at all? One of possible motivations is the desire to stay within the Bayesian framework, which allows us to use integration instead of optimization (a purely technical reason) or assume randomness of parameters (a philosophical reason). This motivation is especially important in complex multiparameter studies, when prior information can be elicited for some of the parameters, but not for the others. If we could express the idea of no prior information as a special noninformative prior distribution, it would allow us to maintain Bayesian setting and derive joint posteriors for our multivariate parameters. For instance, we will be able to use nifty Bayesian software without the necessity to laboriously develop the prior distribution and to justify our subjective point. Just fill the blanks with objective noninformative priors and let computer do the work.

Objective priors are the priors that every expert should agree to use. They are usually weak, nonintrusive, nonviolent priors, which express as little reliance on the prior information as possible and put the emphasis on the data instead. There exist many approaches to the construction of noninformative priors based on certain mathematical principles: Jeffreys prior defined through Fisher information is justified by the idea of invariance to reparameterization [13]; reference priors introduced by Berger and Bernardo maximize the impact of the data on the posterior [2, 4]; maximum entropy priors introduced by Jaynes [12] and Zellner are derived from information-theoretical considerations; and probability matching priors use the idea that Bayesian credible intervals in the absence of prior information should be closely related to the classical confidence intervals.

Bayesian procedures with objective noninformative priors derived by all of these methods bring about results in most cases consistent with the procedures of classical (frequentist) statistics. However the choice of objective prior is not always obvious and sometimes no agreement is achieved between different constructions. Let us use the illustration from the previous subsection.

Ten Coin Tosses

Let us say the coin expert has no past experience with similar coins, and therefore has no opinion at all regarding the probability of heads θ prior to the experiment. Then the envelope is pushed to the statistician (that is you) to decide what prior to put in. An alternative would be to forget about Bayesian approach and resort to MLE, which for 7 observed heads out of 10 tosses provided us in Section 2.4 with the intuitively plausible estimate . However, if we still want to follow the Bayesian approach and take the posterior mean as the estimate of θ, we need to specify the prior, which in our case should be noninformative. An intuitively attractive uniform prior Beta(1, 1) happens to be off. Jeffreys’ prior, reference prior approach, and probability matching prior all point at a different choice of Beta(0.5, 0.5) with a horseshoe-shaped p.d.f

This choice of the horseshoe distribution, though a bit strange at first sight, is justifiable not just mathematically (three different approaches lead to the same decision!) but also intuitively. Just think outside of the coin framework. In general, we are talking about some random event of which we know nothing in advance. What can we say about the probability of this event? Should it be uniformly spread on the interval from 0 to 1? Is it equally likely to be almost impossible (probability close to 0), 50-50 like a fair coin toss (close to the middle of the interval), or almost inevitable (probability close to 1?) Where do the majority of random events in our lives belong?

One can argue that most of the random events in our everyday lives are either almost impossible or almost inevitable. The authors of this book, as most of the other people, have daily routines they tend to follow. Coffee in the morning? Almost inevitable for one of us and almost impossible for the other. Walking to work and back? Almost impossible for one of us and almost inevitable for the other. Both of us show up for classes according to our schedules. There are always slight chances for the opposite to happen, but most of the daily routines are set up firmly. Sun rises, traffic flows, birds fly south in the Fall. You do not often toss a coin to determine whether to eat a lunch or whether to teach a class. If we look closely, there are many random events in our everyday lives, but very few of them have 50-50 chances to happen. Visualize a horseshoe-shaped p.d.f. with the minimum in the middle and most of the probability mass concentrated at both ends.

If you are satisfied with the justification above, let us adopt Beta(0.5, 0.5) prior and evaluate the posterior as we did for subjective prior:

The numerical value of the objective Bayes estimate is

(2.38)

If you still insist on the uniform (you could), this is still a possibility, even if it lacks a good mathematical justification. If taken as a prior, it will not lead us to a grave mistake. We obtain the posterior

and

(2.39)

You can see both priors (uniform and horseshoe) and both corresponding posteriors in Figure 2.11. You can see that the posteriors are relatively close. It can be easily seen that for larger dataset sizes the role of prior becomes less and less visible, and all posteriors look alike. By the way, it is also true for subjective priors.

As we see from Figure 2.11, both Beta(0.5, 0.5) and uniform Beta(1, 1) are not completely uninformative: objective Bayes estimates are still numerically different from MLE. In order to assign a prior for which posterior mean is exactly equal to MLE, we will have to consider an even more radical option Beta(0, 0), also known as Haldane’s prior. This prior is improper. It means that it is not a proper probability distribution in the sense that the area under its density does not integrate to 1 and it even admits events of infinite probability. Though it seems a horrible thing to happen, improper priors play their role in Bayesian analysis, and their use does not bring about any big problems unless they lead to improper posteriors, which is a real disaster preventing any meaningful Bayesian inference. Haldane’s prior with nonempty dataset n ≥ 1 leads to proper posteriors, thus it is usable. We will not emphasize the use of improper priors though sometimes they provide valid noninformative choices.

2.6.3 Empirical Bayes

The empirical Bayes approach to the development of priors is all about using empirical information. In general, we should not use the information from our data to choose the prior. That would constitute a clear case of double dipping: One should not use the same data twice—to choose a prior first, and then to develop the posterior. However, if additional information regarding the distribution of parameter is available, it can be used even if it is not directly relevant to our experiment [9, 19].

Ten Coin Tosses

Suppose we have never tossed our coin yet, so the experiment was not yet performed. It is possible though to have a record of five similar experiments with five other coins which have been tossed before. For the sake of simplicity, assume that each of these five coins was tossed 10 times leading to different numbers of heads: 5, 3, 6, 7, 5. These coins are different from ours, so these results are not directly relevant to the value of θ we want to estimate. However, if we believe the five coins in prior experiments to be somewhat similar to ours, they give us an idea of the distribution of θ over the population of all coins.

If we still decide to go for a beta prior, we will use neither expert estimates nor objective priors, but will estimate prior parameters α and β directly from the empirical data of five prior experiments. We can obtain method of moments estimates

and then using (1.26) solve for and :

Then the empirical Bayes prior is θ ∼ Beta(5.38, 4.97) in Figure 2.12.

To obtain the Bayes estimate for θ for the experiment with 7 heads in 10 tosses, , calculate the posterior

and the estimate for θ is given by

(2.40)

The empirical Bayes approach is more powerful than this example may suggest. First of all, in the presence of empirical information it is not always necessary even to define the prior distribution as belonging to a certain parametric family. Instead, one can use nonparametric Bayes approach.

Nonparametric Bayes modeling is a very hot topic in modern statistics. A Bayesian nonparametric model is by definition a Bayesian model on an infinite-dimensional parametric space, which represents the set of all possible solutions for a given problem (such as the set of all densities in nonparametric density estimation). Bayesian nonparametric models have recently been applied to a variety of machine learning problems including regression, classification, and sequential modeling. However this is out of the focus of the book. We will illustrate nonparametric empirical Bayes on the classical compound Poisson example by Robbins [18], see also Carlin and Louis [8].

Compound Poisson Sample

Let us consider an integer-valued sample x = (x₁, …, x_n) from compound Poisson distribution, where each element has Poisson distribution with its own individual parameter value θ_i corresponding to the distribution mean. This sample can represent a heterogeneous population of insurance customers with Poisson claims rates that vary across the population. Our goal is to estimate the entire vector of parameters θ = (θ₁, …, θ_n). Bayesian approach suggests calculating posterior means

(2.41)

and we will try to do it without explicitly specifying the prior π(θ). Recall that the integral in the denominator of (2.41) is known as the marginal distribution of the data point x_i:

and manipulating the formula of Poisson likelihood

we can also rewrite the numerator in (2.41) as

Combining expressions for the numerator and the denominator, obtain

The only way the prior is included in this formula is inside the marginals. But the marginals can be estimated directly from the data. Suppose that the initial sample contains only k ≤ n distinct values y_j, so that x = (y₁, …, y_k), is the number of x_i: x_i = y_j, and n = ∑^k_{j = 1}N(y_j). Then it is reasonable to use

to estimate the Poisson rates of individual observations.

Nonparametric Bayes approach is becoming increasingly popular in the recent years, especially in the context of discrete distributions for relatively small samples with a large number of parameters.

2.7.1 Exponential Family

The condition of existence of a conjugate prior (2.42) is satisfied for the exponential family, see [7]. Distribution density function for the exponential family can be represented as

(2.42)

Many well-known distributions belong to the exponential family (do not confuse with the exponential distribution defined in Section 1.4: the latter is just a particular case of the exponential family).

The likelihood has a similar form

(2.43)

but the data are fixed and the parameter is treated as a random variable. Thus the factor ∏ⁿ_{i = 1}B(x_i) is a constant with no further effect. The conjugate prior has the same form as the likelihood:

(2.44)

Here k and l are some constants determining the shape of the distribution. Since posterior is proportional to the product of the prior and likelihood, we obtain

(2.45)

Therefore, the posterior belongs to the same conjugate family as the prior with constants k* = k + n; l* = l + ∑ⁿ_{i = 1}T(x_i) depending on the likelihood. Let us provide one example of a conjugate family.

2.7.2 Poisson Likelihood

Poisson distribution introduced in Section 1.2 characterizes random number X of events in a stationary flow occurring during a unit of time. Let us denote the intensity parameter (average number of events during the unit of time) by θ. Distribution of this random variable can be presented as

Poisson distribution belongs to the exponential family with

From this expression, the conjugate prior distribution is

which corresponds to Gamma(α, λ) family with k = λ, l = α − 1.

Let us arrange a sample of n independent copies of the variable X and denote it by (x₁, …, x_n). Likelihood function then can be expressed as

Therefore the posterior can be represented as

which is also a gamma distribution with parameters α + ∑x_i and λ + n. Using formulas for the moments of gamma distribution (1.19), we also obtain the general form of Bayes estimator as the posterior mean:

The smooth transition from prior to posterior in the problem of ten coin tosses discussed in Sections 2.3, 2.4, and 2.6 is explained by conjugate beta prior for binomial likelihood.

2.7.3 Table of Conjugate Distributions

Let us list the basic conjugate families of distributions.