on the probabilistic and statistical apparatus underlying the GUM and related documents

November 9, 2011

1 Preamble

2 Probability

2.1 Meaning

2.2 Chance and Propensity

2.3 Credence and Belief

2.4 Epistemic Probability

2.5 Difficulties

2.6 Role of Context

2.7 example: Prospecting

3 Probability Calculus

3.1 Axioms

3.2 Independence

3.3 Extending the Conversation

3.4 Coherence

3.5 Bayes’s Formula

3.6 example: Tuberculin Test

3.7 Growth of Knowledge

4 Random Variables and Probability Distributions

4.1 Random Variables

4.2 example: Roulette

4.3 example: Light Bulb

4.4 Notational Convention

4.5 Probability Distributions

4.6 Probability Distribution Function

4.7 Probability Density Function

4.8 Expected Value, Variance, and Standard Deviation

4.9 example: Poisson Distribution

4.10 example: Exponential Distribution

4.11 Joint, Marginal, and Conditional Distributions

4.12 example: Shark’s Fin

4.13 Independent Random Variables

4.14 example: Unit Circle

4.15 Correlations and Copulas

5 Functions of Random Variables

5.1 Overview

5.2 Delta Method

5.3 Delta Method — Degeneracy

5.4 example: Radiant Power

5.5 Delta Method — Multivariate

5.6 example: Beer-Lambert-Bouguer Law

5.7 example: Darcy’s Law

5.8 Monte Carlo Method

5.9 example: Volume of Cylinder

5.10 Change-of-Variable Formula — Univariate

5.11 example: Oscillating Mirror

5.12 Change-of-Variable Formula — Multivariate

5.13 example: Linear Combinations of Gaussian Random Variables

5.14 example: Ratio of Exponential Lifetimes

5.15 example: Polar Coordinates

6 Statistical Inference

6.1 Bayesian Inference

6.2 Prior Distribution

6.3 Likelihood Function

6.4 Posterior Distribution

6.5 example: Viral Load

6.6 example: Sleep Hours

6.7 example: Hurricanes

6.8 example: Hypothyroidism

7 Acknowledgments

2 Probability

2.1 Meaning

2.2 Chance and Propensity

2.3 Credence and Belief

2.4 Epistemic Probability

2.5 Difficulties

2.6 Role of Context

2.7 example: Prospecting

3 Probability Calculus

3.1 Axioms

3.2 Independence

3.3 Extending the Conversation

3.4 Coherence

3.5 Bayes’s Formula

3.6 example: Tuberculin Test

3.7 Growth of Knowledge

4 Random Variables and Probability Distributions

4.1 Random Variables

4.2 example: Roulette

4.3 example: Light Bulb

4.4 Notational Convention

4.5 Probability Distributions

4.6 Probability Distribution Function

4.7 Probability Density Function

4.8 Expected Value, Variance, and Standard Deviation

4.9 example: Poisson Distribution

4.10 example: Exponential Distribution

4.11 Joint, Marginal, and Conditional Distributions

4.12 example: Shark’s Fin

4.13 Independent Random Variables

4.14 example: Unit Circle

4.15 Correlations and Copulas

5 Functions of Random Variables

5.1 Overview

5.2 Delta Method

5.3 Delta Method — Degeneracy

5.4 example: Radiant Power

5.5 Delta Method — Multivariate

5.6 example: Beer-Lambert-Bouguer Law

5.7 example: Darcy’s Law

5.8 Monte Carlo Method

5.9 example: Volume of Cylinder

5.10 Change-of-Variable Formula — Univariate

5.11 example: Oscillating Mirror

5.12 Change-of-Variable Formula — Multivariate

5.13 example: Linear Combinations of Gaussian Random Variables

5.14 example: Ratio of Exponential Lifetimes

5.15 example: Polar Coordinates

6 Statistical Inference

6.1 Bayesian Inference

6.2 Prior Distribution

6.3 Likelihood Function

6.4 Posterior Distribution

6.5 example: Viral Load

6.6 example: Sleep Hours

6.7 example: Hurricanes

6.8 example: Hypothyroidism

7 Acknowledgments

Familiarity with the basic concepts and techniques from probability theory and mathematical statistics can best be gained by studying suitable textbooks and exercising those concepts and techniques by solving instructive problems. The books by DeGroot and Schervish [2011], Hoel et al. [1971a], Hoel et al. [1971b], Feller [1968], Lindley [1965a], and Lindley [1965b] are highly recommended for this purpose and appropriate for readers who will have studied mathematical calculus in university courses for science and engineering majors.

This document aims to provide an overview of some of these concepts and techniques that have proven useful in applications to the characterization, propagation, and interpretation of measurement uncertainty as described, for example, by Morgan and Henrion [1992] and Taylor and Kuyatt [1994], and in guidance documents produced by international organizations, including the Guide to the expression of uncertainty in measurement (GUM) [Joint Committee for Guides in Metrology, 2008a] and its supplements [Joint Committee for Guides in Metrology, 2008b]. However, the examples do not all necessarily show a direct connection to measurement science.

Our basic premises are that probability is best suited to express uncertainty quantitatively, and that Bayesian statistical methods afford the best means to exploit information about quantities of interest that originates from multiple sources, including empirical data gathered for the purpose, and preexisting expert knowledge.

Although there is nothing particularly controversial about the calculus of probability or about the mathematical methods of statistics, both the meaning of probability and the interpretation of the products of statistical inference continue to be subjects of debate.

This debate is meta-probabilistic and meta-statistical, in the same sense as metaphysics employs methods different from the methods of physics to study the world. In fact, the debate is liveliest and most scholarly among professional philosophers [Fitelson, 2007]. However, probabilists and statisticians often participate in it when they take off their professional hats and become philosophers [Neyman, 1977], as any inquisitive person is wont to do, at one time or another.

For this reason, we begin with an overview of some of the meanings that have been assigned to probability (§2) before turning to the calculus of probability (§3). In applications, the devices of this calculus are typically brought into play when considering random variables and probability distributions (§4), in particular to characterize the probability distribution of functions of random variables (§5). Statistical inference (§6) uses all of these devices to produce probabilistic statements about unknown quantity values.

In Truth and Probability [Ramsey, 1926, 1931], Frank Ramsey takes the view that probability is “a branch of logic, the logic of partial belief and inconclusive argument”. In this vein, and more generally, probabilities serve to quantify uncertainty. For example, when one states that, with $99\phantom{\rule{0.3em}{0ex}}\%$ confidence, the distance between two geodesic marks is within 0.07 m of 936.84 m, one believes that the actual distance most likely lies between 936.77 m and 936.91 m, but still entertains the possibility that it may lie elsewhere. Similarly, a weather service announcement of $20\phantom{\rule{0.3em}{0ex}}\%$ chance of rain tomorrow for a particular region summarizes an assessment of uncertainty about what will come to pass.

Although relevant to the interpretation of measurement uncertainty, and generally to all applications of probability and statistics, the meaning of probability really is a philosophical issue [Gillies, 2000, Hájek, 2007, Mellor, 2005]. And while there is much disagreement about what probabilities mean, and how they are created to begin with (interpretation and elicitation of probability), there also is essentially universal agreement about how numerical assessments of probability should be manipulated and combined (calculus of probability).

Chances arise in connection with games of chance, and with phenomena that conceivably can recur under essentially the same circumstances. Thus one speaks of the chances of a pair of Kings in a poker hand, or of the chances that the nucleus of an atom of a particular uranium isotope will emit an alpha particle within a given time interval, or of the chances that a person born in France will have blood of type AB. Chances seem to be intrinsic properties of objects or processes in specific environments, maybe propensities for something to happen: their most renowned theorists have been Hans Reichenbach [Reichenbach, 1949], Richard von Mises [von Mises, 1981], and Karl Popper [Popper, 1959].

Credences measure subjective beliefs. They are best illustrated in relation with betting on the outcomes of events one is uncertain about. For example, in this most memorable of bets offered when the thoroughbred Eclipse was about to run against Gower, Tryal, Plume, and Chance in the second heat of the races on May 3rd, 1769, at Epsom Downs: “Eclipse first, the rest nowhere”, with odds of 6-to-4 [Clee, 2007].

The strength or degree of these beliefs can be assessed numerically by techniques that include the observation of betting behavior (actual or virtual), and this assessment can be gauged, and improved, by application of scoring rules [Lindley, 1985].

Dennis Lindley [Lindley, 1985] suggests that degrees of belief can be measured by comparison with a standard, similarly to how length or mass are measured. In general, subjective probabilities can be revealed by judicious application of elicitation methods [Garthwaite et al., 2005].

Consider an urn that contains 100 balls that are identical but for their colors: $\beta $ are black and $100-\beta $ are white. The urn’s contents are thoroughly mixed, and the standard is the probability of the event $\mathit{B}$ of drawing a black ball. Now, given an event $\mathit{E}$, for example that, somewhere in London, it will rain tomorrow, whose probability he wishes to gauge, Peter will select a value for $\beta $ such that he regards gambling on $\mathit{B}$ as equivalent to gambling on $\mathit{E}$ (for the same prize): in these circumstances, $\beta \u2215100$ is Peter’s credence on $\mathit{E}$.

The beliefs that credences measure are subjective and personal, hence the probabilities that gauge them purport to a relationship between a particular knowing subject and the object of this subject’s interest. These beliefs certainly are informed by such knowledge as one may have about a situation, but they also are tempered by one’s preferences or tastes, and do not require that a link be drawn explicitly between that knowledge or sentiment and the corresponding bet. Bruno de Finetti [de Finetti, 1937, 1990], Jimmie Savage [Savage, 1972], and Dennis Lindley [Lindley, 2006] have been leading developers of the subjective, personalistic viewpoint.

Logical (or epistemic, that is, involving or relating to knowledge) probabilities measure the degree to which the truth of a proposition justifies, warrants, or rationally supports the truth of another [Carnap, 1962]. For example, when a medical doctor concludes that a positive result in a tuberculin sensitivity test indicates tuberculosis with $62.5\phantom{\rule{0.3em}{0ex}}\%$ probability (§3.6), or when measurements made during a total eclipse of the sun overwhelmingly favor Einstein’s theory of gravitation over Newton’s [Dyson et al., 1920].

The fact that scientists or judges may not necessarily or explicitly use probabilities to convey their confidence in theories or in arguments [Glymour, 1980] does not reduce the value that probabilities have in models for the rational process of learning from experience, either for human subjects or for reasoning machines that are programmed to make decisions in situations of uncertainty. In this fashion, probability is an extension of deductive logic, and measures degree of confirmation: it does this objectively because it does not involve subjective personal opinion, hence is as incontrovertible as deduction by any of the forms of classical logic.

The difficulty lies in specifying a starting point, a state of a priori ignorance that is similarly objective and hence universally acceptable. Harold Jeffreys [Jeffreys, 1961] provided maybe the first modern, thorough account of how this may be done. He argued, and illustrated in many substantive examples, that it is fit to address the widest range of scientific problems where one wishes to exploit the information in observational data.

The interpretation of probability as an extension of logic makes it particularly well-suited to applications in measurement science, where it is desirable to be able to treat different uncertainty components, which may have been evaluated using different methods, simultaneously, using a uniform vocabulary, and a single set of technical tools. This concept underlies the treatment of measurement uncertainty in the GUM and in its supplements. Richard Cox [Cox, 1946, 1961] and Edwin Jaynes [Jaynes, 1958, 2003] have articulated cogent arguments in support of this view, and José Bernardo [Bernardo, 1979] and James Berger [Berger, 2006] have greatly expanded it.

Even in situations where, on first inspection, chances seem applicable, closer inspection reveals that something else really is needed.

There may be no obvious reason to doubt that the chance is ½ that a coin tossed to assign sides before a football game will land Heads up. However, if the coin instead is spun on its edge on a table, that chance will be closer to either $1\u22153$ or $2\u22153$ than to ½ [Diaconis and Ylvisaker, 1985].

And when it is the magician Persi Warren [DeGroot, 1986] who tosses the coin, then all bets are off because he can manage to toss it so that it always lands Heads up on his hand after he flips it in the air: while the unwary may be willing to bet at even odds on the outcome, for Persi the probability is 1.

And there are situations that naturally lend themselves to, or that even seem to require, multiple interpretations. Take the 20 % chance of rain: does this mean that, of all days when the weather conditions have been similar to today’s in the region the forecast purports to, it has rained some time during the following day with historical frequency of 20 %? Or is this the result of a probabilistic forecast that is to be interpreted epistemically? Maybe it means something else altogether, like: of all things that will fall from the sky tomorrow, 1 in 5 will be a raindrop.

The use of probability as an extension of logic ensures that different people who have the same information (empirical or other) about a measurand should produce the same inferences about it. Example §2.7 illustrates the fact that contextual information relating to a proposition, situation, or event, will influence probabilistic assessments to the extent that different people with different information may, while all acting rationally, produce different uncertainty assessments, and hence different probabilities, for the same proposition, situation, or event.

When probabilities express subjective beliefs, or when they express states of incomplete or imperfect knowledge, different people typically will assign different probabilities to the same statements or events. If they have to reach a consensus on a course of action that is informed by their plural, varied assessments, then they have to engage in a harmonization exercise that preserves the internal coherence of their individual positions. Both statisticians [Stone, 1961, Morris, 1977, Lindley, 1983, Clemen and Winkler, 1999] and philosophers [Bovens and Rabinowicz, 2006, Hartmann and Sprenger, 2011] have addressed this topic.

James and Claire, who both make investments in mining prospects, have been
told that samples from a region surveyed recently have mass fractions of titanium
averaging 3 g kg^{-1}, give or take 1 g kg^{-1} (where “give or take” means that the
true mass fraction of titanium in the region sampled is between 2 g kg^{-1} and
4 g kg^{-1} with 95 % probability). James, however, has also been told that the
samples are of a sandstone with grains of ilmenite. On this basis, James may
assign a much higher probability than Claire to the proposition that asserts that
the region sampled includes an economically viable deposit of titanium
ore.

Once numeric probabilities are in hand, irrespective of how they may be interpreted, the same set of rules, or axioms, is used to combine them. We formulate these axioms in the context where probability is regarded as measuring degree of (rational) belief in the truth of propositions, given a particular body of knowledge and universe of discourse, $\mathit{H}$, that makes all the participating elements meaningful.

Let $\mathit{A}$ and $\mathit{B}$ denote propositions whose probabilities $Pr\left(\mathit{A}|\mathit{H}\right)$ and $Pr\left(\mathit{B}|\mathit{H}\right)$ express degrees of belief about their truth given (or, conditionally upon) the context defined by $\mathit{H}$. The notation $Pr\left(\mathit{B}|\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{H}\right)$ denotes the conditional probability of $\mathit{B}$, assuming that $\mathit{A}$ is true and given the context defined by $\mathit{H}$. Note that $Pr\left(\mathit{B}|\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{H}\right)$ is not necessarily 0 when $Pr\left(\mathit{A}|\mathit{H}\right)=0$: for example, the probability is 0 that a point chosen uniformly at random over the surface of the earth will be on the equator; yet the probability is ½ that, conditionally on its being on the equator, its longitude is between 0° and 180° West of the prime meridian at Greenwich, UK.

The axioms for the calculus of probability are these:

- Convexity:
- $Pr\left(\mathit{A}|\mathit{H}\right)$ is a number between 0 and 1, and it is 1 if and only if $\mathit{H}$ logically implies $\mathit{A}$;
- Addition:
- $Pr\left(\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{or}\phantom{\rule{3.26288pt}{0ex}}\mathit{B}|\mathit{H}\right)=Pr\left(\mathit{A}|\mathit{H}\right)+Pr\left(\mathit{B}|\mathit{H}\right)-Pr\left(\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{B}|\mathit{H}\right)$, where the expression “$\mathit{A}$ or $\mathit{B}$” is true if and only if $\mathit{A}$, $\mathit{B}$, or both are true;
- Multiplication:
- $Pr\left(\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{B}|\mathit{H}\right)=Pr\left(\mathit{B}|\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{H}\right)Pr\left(\mathit{A}|\mathit{H}\right)$.

Since the roles of $\mathit{A}$ and $\mathit{B}$ are interchangeable, the multiplication axiom obviously can also be written as $Pr\left(\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{B}|\mathit{H}\right)=Pr\left(\mathit{A}|\mathit{B}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{H}\right)Pr\left(\mathit{B}|\mathit{H}\right)$.

Most accounts of mathematical probability theory use an additional rule (countable additivity) that ensures that the probability that at least one proposition is true, among countably infinitely many mutually exclusive propositions, equals the sum of their individual probabilities [Casella and Berger, 2002, Definition 1.2.4]. (“Countably infinitely many” means “as many as there are integer numbers”.)

When the context that $\mathit{H}$ defines is obvious, often one suppresses explicit reference to it, as in this derivation: if $\tilde{\mathit{A}}$ denotes the negation of $\mathit{A}$, then Convexity and the Addition Rule imply that $1$ $=Pr\left(\mathit{A}\phantom{\rule{3.26288pt}{0ex}}\text{or}\phantom{\rule{3.26288pt}{0ex}}\tilde{\mathit{A}}\right)$ $=Pr\left(\mathit{A}\right)$ $+Pr\left(\tilde{\mathit{A}}\right)$ because one but not both of $\mathit{A}$ and $\tilde{\mathit{A}}$ must be true.

The concept of independence pervades much of probability theory. Two propositions $\mathit{A}$ and $\mathit{B}$ are independent if the probability that both are true equals the product of their individual probabilities of being true. If $\mathit{A}$ asserts that there is a Queen in Alexandra’s poker hand, and $\mathit{B}$ asserts that Beatrice’s comprises red cards only, both hands having been dealt from the same deck, then $\mathit{A}$ and $\mathit{B}$ are independent. Intuitively, if knowledge of the truth of one proposition influences the assessment of probability of another, then they are dependent: in particular, two mutually exclusive propositions are dependent. If $Pr\left(\mathit{B}|\mathit{A}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\mathit{H}\right)=Pr\left(\mathit{B}|\mathit{H}\right)$, then $\mathit{A}$ and $\mathit{B}$ are independent given $\mathit{H}$.

When considering the probability of a proposition, it often proves advantageous to consider the truth or falsity of another one, somehow related to the first [Lindley, 2006, §5.6]. To assess the probability $Pr\left(+\right)$ of a positive tuberculin skin test (§3.6), it is convenient to consider how the test performs separately in persons infected or not infected with Mycobacterium tuberculosis: if $\mathit{I}$ denotes infection, then $Pr\left(+\right)$ $=Pr\left(+|\mathit{I}\right)Pr\left(\mathit{I}\right)$ $+Pr\left(+|\tilde{\mathit{I}}\right)Pr\left(\tilde{\mathit{I}}\right)$, where $Pr\left(+|\tilde{\mathit{I}}\right)$ is the probability of a false positive, and $1-Pr\left(+|\mathit{I}\right)$ is the probability of a false negative, both more accessible than $Pr\left(+\right)$.

If an ideal reasoning agent (human or machine) assigns probabilities to events or to the truth of propositions according to the foregoing axioms, then this agent’s beliefs are said to be coherent. In these circumstances, if probabilities are used to inform bets concerning the truth of propositions in the universe of discourse where these probabilities are meaningful, then it is impossible (for a “bookie”) to devise a collection of bets that bring an assured loss to this agent (a so-called “Dutch Book”).

Now, suppose that, having ascertained the truth of a proposition $\mathit{A}$, one produces $Pr\left(\mathit{C}|\mathit{A}\right)$ as assessment of $\mathit{C}$’s truth on the evidence provided by $\mathit{A}$. Next, one determines that $\mathit{B}$, too, is true and revises this last assessment of $\mathit{C}$’s truth to become $Pr\left(\mathit{C}|\mathit{B}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{A}\right)$. The process whereby probabilities are updated is coherently extensible if the resulting assessment is the same irrespective of whether the evidence provided by $\mathit{A}$ and $\mathit{B}$ is brought to bear either sequentially, as just considered, or simultaneously. The incorporation of information from multiple sources, and the corresponding propagation of uncertainty, that is carried out by application of Bayes’ formula, which is described next and illustrated in examples §3.6 and §6.7, is coherently extensible.

If exactly one (that is, one and one only) among propositions ${\mathit{A}}_{1},\dots ,{\mathit{A}}_{\mathit{n}}$ can be true, and $\mathit{B}$ is another proposition with positive probability, then

This follows from the axioms above because $Pr\left({\mathit{A}}_{\mathit{j}}|\mathit{B}\right)=Pr\left({\mathit{A}}_{\mathit{j}}\phantom{\rule{3.26288pt}{0ex}}\text{and}\phantom{\rule{3.26288pt}{0ex}}\mathit{B}\right)\u2215Pr\left(\mathit{B}\right)$ (Multiplication axiom), whose numerator equals $Pr\left(\mathit{B}|{\mathit{A}}_{\mathit{j}}\right)Pr\left({\mathit{A}}_{\mathit{j}}\right)$ (Multiplication axiom), and whose denominator equals $Pr\left(\mathit{B}|{\mathit{A}}_{1}\right)Pr\left({\mathit{A}}_{1}\right)$ $\dots \phantom{\rule{0.3em}{0ex}}$ $Pr\left(\mathit{B}|{\mathit{A}}_{\mathit{n}}\right)Pr\left({\mathit{A}}_{\mathit{n}}\right)$ (“extending the conversation”, as in §3.3).

Richard has been advised that his tuberculin skin test has returned a positive result. The tuberculin skin test has a reported false-negative rate of $25\phantom{\rule{0.3em}{0ex}}\%$ during the initial evaluation of persons with active tuberculosis [American Thoracic Society, 1999, Holden et al., 1971]: this means that the probability is 0.25 that the test will yield a negative ($-$) response when administered to an infected person ($\mathit{I}$), $Pr\left(-|\mathit{I}\right)=0.25$. Therefore, the probability is only 0.75 that infection will yield a positive test result. In populations where cross-reactivity with other mycobacteria is common, the test’s false-positive rate is 5%: that is, the conditional probability of a positive result ($+$) for a person that is not infected ($\tilde{\mathit{I}}$) is $Pr\left(+|\tilde{\mathit{I}}\right)=0.05$.

Richard happens to live in an area where tuberculosis has a prevalence of $10\phantom{\rule{0.3em}{0ex}}\%$. Given the positive result of the test he underwent, the probability that he is infected is

$$\begin{array}{llll}\hfill Pr\left(\mathit{I}|+\right)=& \frac{Pr\left(+|\mathit{I}\right)Pr\left(\mathit{I}\right)}{Pr\left(+|\mathit{I}\right)Pr\left(\mathit{I}\right)+Pr\left(+|\tilde{\mathit{I}}\right)Pr\left(\tilde{\mathit{I}}\right)}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \frac{\left(0.75\times 0.10\right)}{\left(0.75\times 0.10\right)+\left(0.05\times 0.90\right)}=0.625.\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$$Common sense suggests that the diagnostic value of the test should depend on its false-negative and false-positive rates, as well as on the prevalence of the disease: Bayes’ formula states exactly how these ingredients should be combined to produce $Pr\left(\mathit{I}|+\right)$, which expresses that diagnostic value quantitatively.

Richard has the tuberculin skin test repeated, and this second test also turns out positive. To incorporate this additional piece of evidence into the probability that Richard is infected, first we summarize the state of knowledge (about whether he is infected) determined by the result from the first test. This is done by defining $\mathit{Q}\left(\mathit{I}\right)=Pr\left(\mathit{I}|+\right)=0.625$ and $\mathit{Q}\left(\tilde{\mathit{I}}\right)=1-\mathit{Q}\left(\mathit{I}\right)=0.375$, and using them in the role that the overall probability of infection ($10\phantom{\rule{0.3em}{0ex}}\%$) or non-infection ($90\phantom{\rule{0.3em}{0ex}}\%$) played prior to Richard’s first test, when all one knew about his condition was that he was a member of a population where the prevalence of tuberculosis was $10\phantom{\rule{0.3em}{0ex}}\%$.

Again applying Bayes’ theorem, and assuming that the two tests are independent, the revised probability that Richard is infected after two positive tests is

$$\begin{array}{llll}\hfill \mathit{Q}\left(\mathit{I}|+\right)=& \frac{Pr\left(+|\mathit{I}\right)\mathit{Q}\left(\mathit{I}\right)}{Pr\left(+|\mathit{I}\right)\mathit{Q}\left(\mathit{I}\right)+Pr\left(+|\tilde{\mathit{I}}\right)\mathit{Q}\left(\tilde{\mathit{I}}\right)}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \frac{\left(0.75\times 0.625\right)}{\left(0.75\times 0.625\right)+\left(0.05\times 0.375\right)}=0.962.\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$$If, instead, one had been initially told that Richard had had two independent, positive tuberculin skin tests, then the calculation would have been:

$$\begin{array}{llll}\hfill Pr\left(\mathit{I}|++\right)=& \frac{Pr\left(++|\mathit{I}\right)Pr\left(\mathit{I}\right)}{Pr\left(++|\mathit{I}\right)Pr\left(\mathit{I}\right)+Pr\left(++|\tilde{\mathit{I}}\right)Pr\left(\tilde{\mathit{I}}\right)}\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill =& \frac{\left(0.7{5}^{2}\times 0.10\right)}{\left(0.7{5}^{2}\times 0.10\right)+\left(0.0{5}^{2}\times 0.90\right)}=0.962.\phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\end{array}$$This example illustrates the fact that Bayes’ theorem produces the same probability irrespective of whether the information is incorporated sequentially, or all at once.

The example in §3.6 illustrated how the probability of Richard being infected increased (relative to the overall probability of infection in the town where he lives) as a first, and then a second tuberculin test turned out positive. However, even if he is infected, by chance alone a test may turn out negative. In a sequence of tests, therefore, the probability of his being infected may oscillate, increasing when a test turns out positive, decreasing when some subsequent test turns out negative.

Therefore, the question naturally arises of whether a person employing the Bayesian method of exploiting information, and incorporating it into the current state of knowledge, ever will, in situations of uncertainty, arrive at conclusions with overwhelming confidence. Jimmie Savage proved rigorously that the answer is “yes” with great generality: “with observation of an abundance of relevant data, the person is almost certain to become highly convinced of the truth, and […] he himself knows this to be the case” [Savage, 1972, §3.6].

The restriction to “relevant data” is critical: in relation with the tuberculin test, if it happened that $Pr\left(+|\mathit{I}\right)=Pr\left(+|\tilde{\mathit{I}}\right)=\text{\xbd}$, then the test would have no discriminatory power, and in fact would be irrelevant to learning about disease status.

The notion of random variable originates in games of chance, like roulette, whose outcomes are unpredictable. Its rigorous mathematical definition (measurable function from one probability space into another) is unlikely to be of great interest to the metrologist. Instead, one may like to keep in mind its heuristic meaning: the value of a quantity that has a probability distribution as an attribute whose role is to describe the uncertainty associated with that value.

In the version of roulette played in Monte Carlo, the possible outcomes are numbers in the set $\mathcal{\mathcal{X}}=\left\{0,1,\dots ,36\right\}$ (usually one disregards other possible, but “uninteresting” outcomes, including those where the ball exits the wheel and lands elsewhere, or where it lands inside the wheel but in none of its numbered pockets). Once those 37 numbers are deemed to be equally likely, one can speak of a random variable that is equal to $0$ with probability $1\u221537$, or that is odd with probability $18\u221537$. (Note that these statements are meaningful irrespective of whether the event in question will happen in the future, or has happened already, provided one does not know its actual outcome yet).

The GE A19 Party Light 25 W incandescent light bulb has expected lifetime 2000 h: this is usually taken to mean that, if a brand new bulb is turned on and left on supplied with constant 120 V electrical current until it burns out, its actual lifetime may be described as a realized value (realization, or outcome) of a random variable with an exponential probability distribution (§4.10) whose expected value is $\text{2000}\phantom{\rule{0.3em}{0ex}}\text{h}$ — this is denoted $\eta $ in §4.10, and in general it needs to be estimated from experimental data.

The concept of random variable applies just as well to domains of discourse unrelated to games of chance, hence can be used to suggest uncertainty about the value of a quantity, irrespective of the source of this uncertainty, including situations where there is nothing “random” (in the sense of “chancy”) in play.

For the most part, upper case letters (Roman or Greek) denote generic quantity values modeled as random variables, and their lowercase counterparts denote particular values.

Upper case letters like $\mathit{X}$ or ${\mathit{X}}_{1},{\mathit{X}}_{2},\dots \phantom{\rule{0.3em}{0ex}}$, and $\mathit{Y}$, denote generic random variables, without implying that any of the former necessarily play the role of input quantity values (as defined in the international vocabulary of metrology (VIM) Joint Committee for Guides in Metrology [2008c], VIM 2.50), or that the latter necessarily plays the role of output quantity values (VIM 2.51) in a measurement model (VIM 2.48).

The probabilities most commonly encountered in metrological practice concern sets of numbers that a quantity value may take: in this case, if $\mathit{X}$ denotes a random variable whose values belong to a set $\mathcal{\mathcal{X}}$, and $\mathit{A}$ is a subset of $\mathcal{\mathcal{X}}$, then $Pr\left(\mathit{X}\in \mathit{A}\right)$ denotes the probability that $\mathit{X}$’s value lies in $\mathit{A}$. For example, if $\mathit{X}$ represents the length (expressed in meter, say) of a gauge block, then $\mathcal{\mathcal{X}}$ would be the set of all possible values of length, and $\mathit{A}$ could be the subset of such values between 0.0423 m and 0.0427 m, say.

Given a random variable $\mathit{X}$ one can then define a function ${\mathit{P}}_{\mathit{X}}$ such that ${\mathit{P}}_{\mathit{X}}\left(\mathit{A}\right)=Pr\left(\mathit{X}\in \mathit{A}\right)$ for all $\mathit{A}\subset \mathcal{\mathcal{X}}$ to which a probability can be assigned. This ${\mathit{P}}_{\mathit{X}}$ is called $\mathit{X}$’s probability distribution.

If $\mathcal{\mathcal{X}}$ is countable (that is, either finite or infinite but with as many elements as there are positive integers), then one says that $\mathit{X}$ has a discrete distribution, which is fully specified by the probability it assigns to each value in $\mathcal{\mathcal{X}}$. For example, the outcome of a roulette wheel is a random variable whose probability distribution is discrete.

If $\mathcal{\mathcal{X}}$ is uncountable (that is, it has as many elements as there are real numbers) and $Pr\left(\mathit{X}=\mathit{x}\right)=0$ for all $\mathit{x}\in \mathcal{\mathcal{X}}$, then one says that $\mathit{X}$ has a continuous distribution. For example, the lifetime of an incandescent light bulb that does light up and then is constantly left on until it burns out is a random variable with a continuous distribution.

A distribution may be neither discrete nor continuous, but of a mixed type instead: for example, when a random variable is equal to 0 with probability $\u03f5>0$, and has an exponential distribution (see §4.10) with probability $1-\u03f5$. Since a brand new light bulb has a positive probability of burning out the instant it is turned on, its lifetime may more realistically be modeled as a random variable that has an “atom” of probability at 0, and is exponential with the complementary probability.

The probability distribution of a random variable $\mathit{X}$ whose possible values are real numbers, can be succinctly described by its probability distribution function, which is the function ${\mathit{P}}_{\mathit{X}}$ such that ${\mathit{P}}_{\mathit{X}}\left(\mathit{x}\right)=Pr\left(\mathit{X}\le \mathit{x}\right)$ for every real number $\mathit{x}$.

Note that the symbol we use here to denote the probability distribution function, is the same that we used in §4.5 to denote the probability distribution itself. Any confusion this may cause will be promptly resolved by examining the argument of ${\mathit{P}}_{\mathit{X}}$: if it is a set, then we mean the distribution itself, while if it is a number or a vector with numerical components, then we mean the probability distribution function.

For example, if $\mathit{X}$ is real-valued and $\mathit{x}$ is a particular real number, then in ${\mathit{P}}_{\mathit{X}}\left(\mathit{x}\right)={\mathit{P}}_{\mathit{X}}\left(\left(-\infty ,\mathit{x}\right]\right)$ the ${\mathit{P}}_{\mathit{X}}$ on the left hand side refers to the probability distribution function, while the ${\mathit{P}}_{\mathit{X}}$ on the right hand side refers to the distribution itself because $\left(-\infty ,\mathit{x}\right]$ denotes the set of all real numbers no greater than $\mathit{x}$. Since the distribution function determines the distribution, the confusion is harmless.

If $\mathit{X}$ has a discrete distribution (§4.5), then its probability density (also known as its probability mass function) is the function ${\mathit{p}}_{\mathit{X}}$ such that ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)=Pr\left(\mathit{X}=\mathit{x}\right)$ for $\mathit{x}\in \mathcal{\mathcal{X}}$.

If $\mathcal{\mathcal{X}}$ is uncountable and $\mathit{X}$’s distribution is continuous and sufficiently smooth (in the sense described next), then the corresponding probability density function (PDF) is defined similarly to a material object’s mass density, as follows.

Consider the simplest case, where $\mathcal{\mathcal{X}}$ is an interval of real numbers, and suppose that $\mathit{x}$ is one point in the interior of this interval. Now suppose that ${\delta}_{1}>{\delta}_{2}>\dots \phantom{\rule{0.3em}{0ex}}$ is an infinite sequence of positive numbers decreasing to zero. If $\mathit{X}$’s probability distribution is sufficiently smooth, then the limit ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)=\underset{\mathit{n}\to \infty}{lim}\left(\right.{\mathit{P}}_{\mathit{X}}\left(\mathit{x}+{\delta}_{\mathit{n}}\right)-{\mathit{P}}_{\mathit{X}}\left(\mathit{x}-{\delta}_{\mathit{n}}\right)\left)\right.\u2215\left(2{\delta}_{\mathit{n}}\right)$ exists. The function ${\mathit{p}}_{\mathit{X}}$ so defined is $\mathit{X}$’s probability density function. If the distribution function is differentiable, then the probability density is the derivative of the probability distribution function, ${\mathit{p}}_{\mathit{X}}={\mathit{P}}_{\mathit{X}}^{\prime}$.

Both the probability distribution function and the probability density function have multivariate counterparts.

The expectation (expected value, or mean value) of a (scalar or vector valued) function $\phi $ of a random variable $\mathit{X}$ is $\mathbb{\mathbb{E}}\left(\phi \left(\mathit{X}\right)\right)={\int}_{\mathcal{\mathcal{X}}}\phi \left(\mathit{x}\right){\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)\text{d}\mathit{x}$ if $\mathit{X}$ has a continuous probability distribution with density ${\mathit{p}}_{\mathit{X}}$, or $\mathbb{\mathbb{E}}\left(\phi \left(\mathit{X}\right)\right)={\sum}_{\mathit{x}\in \mathcal{\mathcal{X}}}\phi \left(\mathit{x}\right){\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)$ if $\mathit{X}$ has a discrete distribution. Note that $\mathbb{\mathbb{E}}\left(\phi \left(\mathit{X}\right)\right)$ can be computed without determining the probability distribution of the random variable $\phi \left(\mathit{X}\right)$ explicitly.

$\mathbb{\mathbb{E}}\left(\mathit{X}\right)$ indicates $\mathit{X}$’s location, or the center of its probability distribution: therefore it is a most succinct summary of this distribution, and it is the best estimate of $\mathit{X}$’s value in the sense that it has the smallest mean squared error.

The median is another indication of location for a scalar random variable: it is any value $\xi $ such that $Pr\left(\mathit{X}\le \xi \right)\ge \text{\xbd}$ and $Pr\left(\mathit{X}\ge \xi \right)\ge \text{\xbd}$, and it need not be unique. The median is the best estimate of $\mathit{X}$’s value in the sense that it has the smallest absolute deviation.

Neither the mean nor the median need be “representative” values of the distribution. For example, when $\mathit{X}$ denotes a proportion whose most common values are close to 0 or to 1 and its mean is close to ½, then values close to the mean are very unlikely. $\mathbb{\mathbb{E}}\left(\mathit{X}\right)$ need not exist (in the sense that the defining integral or sum may fail to converge).

$\mathbb{\mathbb{E}}\left({\mathit{X}}^{\mathit{k}}\right)$, where $\mathit{k}$ is a positive integer, is called $\mathit{X}$’s $\mathit{k}$th moment. The variance of $\mathit{X}$ is ${\sigma}^{2}=\mathbb{\mathbb{V}}\left(\mathit{X}\right)=\mathbb{\mathbb{E}}\left[{\left(\mathit{X}-\mathbb{\mathbb{E}}\left(\mathit{X}\right)\right)}^{2}\right]$, or, equivalently, the difference between its second moment and the square of its first moment. The positive square root of the variance, $\sigma $, is the standard deviation of $\mathit{X}$.

The only values that a Poisson distributed random variable $\mathit{X}$ can take are the non-negative integers: 0, 1, 2, …, and the probability that its value is $\mathit{x}$ is ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)={\lambda}^{\mathit{x}}{\mathit{e}}^{-\lambda}\u2215\mathit{x}!$, where $\lambda $ is some given positive number, and $\mathit{x}!=\mathit{x}\left(\mathit{x}-1\right)\dots 1$. This model distributes its unit of probability into infinitely many lumps, one at each non-negative integer, so that ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)$ decreases rapidly with increasing $\mathit{x}$, and ${\mathit{p}}_{\mathit{X}}\left(0\right)+{\mathit{p}}_{\mathit{X}}\left(1\right)+{\mathit{p}}_{\mathit{X}}\left(2\right)+\cdots =1$. Both the expected value and the variance equal $\lambda $. The number of alpha particles emitted by a sample containing the radionuclide ${}^{210}\text{Po}$, during a period of $\mathit{t}$ seconds that is a small fraction of this isotope’s half-life (138 days), is a value of a Poisson random variable with mean proportional to $\mathit{t}$.

Suppose that $\mathit{X}$ represents the lifetime (thousands of hours) of an incandescent light bulb, such that, for $0<\mathit{a}<\mathit{b}$, $Pr\left(\mathit{a}<\mathit{X}<\mathit{b}\right)=exp\left(-\mathit{a}\u2215\eta \right)-exp\left(-\mathit{b}\u2215\eta \right)$, for some given number $\eta >0$: note that, as $\mathit{a}$ decreases toward 0, and $\mathit{b}$ increases without limit, $Pr\left(\mathit{a}<\mathit{X}<\mathit{b}\right)$ approaches 1. Focus on a particular number $\mathit{x}>0$, and consider the ratio $Pr\left(\mathit{x}-\delta <\mathit{X}<\mathit{x}+\delta \right)\u2215\left(2\delta \right)=exp\left(-\mathit{x}\u2215\eta \right)\left[exp\left(\delta \u2215\eta \right)-exp\left(-\delta \u2215\eta \right)\right]\u2215\left(2\delta \right)$ for some $\delta >0$. As $\delta $ decreases to $0$ this ratio approaches $\left(1\u2215\eta \right)exp\left(-\mathit{x}\u2215\eta \right)$. Therefore, the function ${\mathit{p}}_{\mathit{X}}$ such that ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)=\left(1\u2215\eta \right)exp\left(-\mathit{x}\u2215\eta \right)$ is the probability density of the exponential distribution. In this case, the probability distribution function is ${\mathit{P}}_{\mathit{X}}$ such that ${\mathit{P}}_{\mathit{X}}\left(\mathit{x}\right)=Pr\left(\mathit{X}\le \mathit{x}\right)=1-exp\left(-\mathit{x}\u2215\eta \right)$. $\mathit{X}$’s mean value is $\mathbb{\mathbb{E}}\left(\mathit{X}\right)=\eta $, and its variance is $\mathbb{\mathbb{V}}\left(\mathit{X}\right)={\eta}^{2}$. Figure 1 illustrates both the distribution function and the density for the case where $\eta =\text{2000}\phantom{\rule{0.3em}{0ex}}\text{h}$.

Suppose that $\mathit{X}$ represents a bivariate quantity value, for example, the Cartesian coordinates $\left(\mathit{U},\mathit{V}\right)$ of a point inside a circle of unit radius centered at $\left(0,0\right)$. In this case the range $\mathcal{\mathcal{X}}$ of $\mathit{X}=\left(\mathit{U},\mathit{V}\right)$ is this unit circle. The joint probability distribution of $\mathit{U}$ and $\mathit{V}$ describes a state of knowledge about the location of $\mathit{X}$: for example, that more likely than not $\mathit{X}$ is less than ½ away from the center of the circle: statements of this kind involve $\mathit{U}$ and $\mathit{V}$ together (that is, jointly).

The marginal probability distributions of $\mathit{U}$ and $\mathit{V}$ are the probability distributions that characterize the state of knowledge about each of them separately from the other: for example, that more likely than not $-\text{\xbd}<\mathit{U}<\text{\xbd}$, irrespective of $\mathit{V}$.

Clearly, the marginal distributions have to be consistent with the joint distribution, and while it is true that the joint distribution determines the marginal distributions, the reverse is not true, in that typically there are many joint distributions consistent with given marginal distributions [Possolo, 2010].

Now, suppose one knows that $\mathit{U}=2\u22153$. This implies that $-\sqrt{5}\u22153<\mathit{V}<\sqrt{5}\u22153$, hence that $\mathit{X}=\left(\mathit{U},\mathit{V}\right)$ is somewhere on a particular chord $\mathcal{\mathcal{C}}$ of the unit circle. The conditional probability distribution of $\mathit{V}$ given that $\mathit{U}=2\u22153$ is a (univariate) probability distribution over this chord.

The random variables $\mathit{X}$ and $\mathit{Y}$ take values in the interval $\left(1,2\right)$ and have joint probability density function ${\mathit{p}}_{\mathit{X},\mathit{Y}}$ such that ${\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)=\left(\mathit{x}+\mathit{y}\right)\u22153$ for $1\le \mathit{x},\mathit{y}\le 2$ and is zero otherwise (Figure 2): since ${\mathit{p}}_{\mathit{X},\mathit{Y}}\ge 0$ and ${\int}_{1}^{2}{\int}_{1}^{2}{\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)\text{d}\mathit{x}\text{d}\mathit{y}=1$, ${\mathit{p}}_{\mathit{X},\mathit{Y}}$ is a bona fide (bivariate) probability density.

The density of the marginal distribution of $\mathit{X}$ is ${\mathit{p}}_{\mathit{x}}\left(\mathit{x}\right)={\int}_{1}^{2}{\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)\text{d}\mathit{y}=\mathit{x}\u22153+1\u22152$, and similarly for $\mathit{Y}$. And the density of the conditional distribution of $\mathit{Y}$ given $\mathit{X}=\mathit{x}$ is ${\mathit{p}}_{\mathit{Y}|\mathit{X}}\left(\mathit{y}|\mathit{x}\right)={\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)\u2215{\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)=\left(\mathit{x}+\mathit{y}\right)\u2215\left(\mathit{x}+3\u22152\right)$.

To determine the probability density of $\mathit{R}=\mathit{Y}\u2215\mathit{X}$, also depicted in Figure 2, note that $\text{\xbd}\le \mathit{R}\le 2$. First, consider the case $\text{\xbd}\mathit{r}\le 1$: $Pr\left(\mathit{R}\le \mathit{r}\right)={\left(2\mathit{r}-1\right)}^{2}\u2215\left(2\mathit{r}\right)$ and ${\mathit{p}}_{\mathit{R}}\left(\mathit{r}\right)=\left(4{\mathit{r}}^{2}-1\right)\u2215\left(2{\mathit{r}}^{2}\right)$. Next, consider the case $1<\mathit{r}\le 2$: $Pr\left(\mathit{R}\le \mathit{r}\right)=1-{\left(2-\mathit{r}\right)}^{2}\u2215\left(2\mathit{r}\right)$ and ${\mathit{p}}_{\mathit{R}}\left(\mathit{r}\right)=\left(4-{\mathit{r}}^{2}\right)\u2215\left(2{\mathit{r}}^{2}\right)$.

The (scalar or vectorial) random variables $\mathit{X}$ and $\mathit{Y}$ are independent if and only if $Pr\left(\mathit{X}\in \mathit{A}\phantom{\rule{1em}{0ex}}\text{and}\phantom{\rule{1em}{0ex}}\mathit{Y}\in \mathit{B}\right)$ $=Pr\left(\mathit{X}\in \mathit{A}\right)Pr\left(\mathit{Y}\in \mathit{B}\right)$ for all subsets $\mathit{A}$ and $\mathit{B}$ in their respective ranges (to which probabilities can be coherently assigned.) Suppose $\mathit{X}$ and $\mathit{Y}$ have joint probability distribution with probability density function ${\mathit{p}}_{\mathit{X},\mathit{Y}}$, and marginal density functions ${\mathit{p}}_{\mathit{X}}$ and ${\mathit{p}}_{\mathit{Y}}$: the random variables are independent if and only if ${\mathit{p}}_{\mathit{X},\mathit{Y}}={\mathit{p}}_{\mathit{X}}{\mathit{p}}_{\mathit{Y}}$.

Suppose that the probability distribution of a point is uniform inside the circle of unit radius centered at the origin $\left(0,0\right)$ of the Euclidean plane. This means that the probability that a point with Cartesian coordinates $\left(\mathit{X},\mathit{Y}\right)$ should lie in a subset $\mathit{S}$ of this circle is proportional to $\mathit{S}$’s area, but is otherwise independent of $\mathit{S}$’s shape or location within the circle. The probability density function of the joint distribution of $\mathit{X}$ and $\mathit{Y}$ is the function ${\mathit{p}}_{\mathit{X},\mathit{Y}}$ such that ${\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)=1\u2215\pi $ if ${\mathit{x}}^{2}+{\mathit{y}}^{2}<1$, and ${\mathit{p}}_{\mathit{X},\mathit{Y}}\left(\mathit{x},\mathit{y}\right)=0$ otherwise. The random variables $\mathit{X}$ and $\mathit{Y}$ are dependent (§4.13): for example, if one is told that $\mathit{X}>\text{\xbd}$, then one can surely conclude that $-\sqrt{3}\u22152<\mathit{Y}<\sqrt{3}\u22152$. The marginal distribution of $\mathit{X}$ has density ${\mathit{p}}_{\mathit{X}}$ such that ${\mathit{p}}_{\mathit{X}}\left(\mathit{x}\right)=\left(2\u2215\pi \right){\int}_{0}^{\sqrt{1-{\mathit{x}}^{2}}}\text{d}\mathit{y}=\left(2\u2215\pi \right)\sqrt{1-{\mathit{x}}^{2}}$ for $-1<\mathit{x}<1$. $\mathit{X}$ has expected value 0 and standard deviation ½. Owing to symmetry, $\mathit{X}$ and $\mathit{Y}$ have identical marginal distributions.

If two or more of the random variables are dependent, then modeling their individual probability distributions will not suffice to specify their joint behavior: their joint probability distribution is needed.

One commonly used metric of dependence between two random variables $\mathit{X}$ and $\mathit{Y}$ is Pearson’s product-moment correlation coefficient, defined as $\rho \left(\mathit{X},\mathit{Y}\right)=\mathbb{\mathbb{E}}\left[\left(\mathit{X}-\mathbb{\mathbb{E}}\left(\mathit{X}\right)\right)\left(\mathit{Y}-\mathbb{\mathbb{E}}\left(\mathit{Y}\right)\right)\right]\u2215\sqrt{\mathbb{\mathbb{V}}\left(\mathit{X}\right)\mathbb{\mathbb{V}}\left(\mathit{Y}\right)}$. However, it is possible for the variables to be dependent and still have $\rho \left(\mathit{X},\mathit{Y}\right)=0$.

When the only information in hand are the expected values, standard deviations, and correlations, and still one needs a joint distribution consistent with this information, then the usual course of action is to assign distributions to them individually, and then manufacturing a joint distribution using a copula [Possolo, 2010] — there is, however, a multitude of different copulas that can be used for this purpose, and the choice that must be made generally is influential.

If a random variable $\mathit{Y}$ is a function of other random variables, $\mathit{Y}=\phi \left({\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}\right)$, then $\phi $ and the joint probability distribution of ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ determine the probability distribution of $\mathit{Y}$.

If only the means, standard deviations, and correlations of ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ are known, then it still is possible to derive approximations to the mean and standard deviation of $\mathit{Y}$, by application of the Delta Method.

If the joint probability distribution of ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ is known, then it may be possible to determine the probability distribution of $\mathit{Y}$ analytically, using the change of variables formula.

In general, it is possible to obtain a sample from $\mathit{Y}$’s distribution by taking a sample from the joint distribution of ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ and applying $\phi $ to each element of this sample (§5.8). The results may then be summarized in several different ways: one of them is an estimate of the probability density of $\mathit{Y}$ [Silverman, 1986], a procedure that is implemented in function density of the R environment for statistical programming and graphics [R Development Core Team, 2010].

If $\mathit{X}$ is a random variable with mean $\mu $ and variance ${\sigma}^{2}$, $\phi $ is a differentiable real-valued function of a real variable whose first derivative does not vanish at $\mu $, and $\mathit{Y}=\phi \left(\mathit{X}\right)$, then $\mathbb{\mathbb{E}}\left(\mathit{Y}\right)\approx \phi \left(\mu \right)$, and $\mathbb{\mathbb{V}}\left(\mathit{Y}\right)\approx {\left[{\phi}^{\prime}\left(\mu \right)\right]}^{2}{\sigma}^{2}$. (This results from the so-called Taylor approximation that replaces $\phi $ by a straight line tangent to its graph at $\mu $.)

If $\mathit{X}=\left({\mathit{V}}_{1}+\cdots +{\mathit{V}}_{\mathit{m}}\right)\u2215\mathit{m}$ is an average of independent, identically distributed random variables with finite variance, then $\sqrt{\mathit{m}}\left(\phi \left({\mathit{V}}_{\mathit{m}}\right)-\phi \left(\mu \right)\right)$ also is approximately Gaussian with mean $0$ and standard deviation $\left|{\phi}^{\prime}\left(\mu \right)\right|\sigma $, where $\left|{\phi}^{\prime}\left(\mu \right)\right|$ denotes the absolute value of the first derivative of $\phi $ evaluated at $\mu $. The quality of the approximation improves with increasing $\mathit{m}$.

When ${\phi}^{\prime}\left(\mu \right)=0$ and ${\phi}^{\prime \prime}\left(\mu \right)$ exists and is not zero, $\mathit{X}=\left({\mathit{V}}_{1}+\cdots +{\mathit{V}}_{\mathit{m}}\right)\u2215\mathit{m}$ is an average of independent, identically distributed random variables with finite variance, and $\mathit{m}$ is large, then the probability distribution of $\mathit{m}\left(\phi \left(\mathit{X}\right)-\phi \left(\mu \right)\right)$ is approximately like that of ${\sigma}^{2}{\phi}^{\prime \prime}\left(\mu \right){\mathit{Z}}^{2}\u22152$, where $\mathit{Z}$ denotes a Gaussian (or, normal) random variable with mean 0 and standard deviation 1. Since the variance of ${\mathit{Z}}^{2}$ is 2, the standard deviation of $\phi \left({\mathit{V}}_{\mathit{m}}\right)$ is approximately ${\sigma}^{2}\left|{\phi}^{\prime \prime}\left(\mu \right)\right|\u2215\sqrt{2}$, rather different from what applies in the conditions of §5.2.

Consider a surface whose reflectance is Lambertian: that is, light falling on it is scattered in such a way that the surface’s brightness apparent to an observer is the same regardless of the observer’s angle of view. The radiant power $\mathit{W}$ emitted by such surface that is measured by a sensor aimed at angle $\mathit{A}$ to the surface’s normal is proportional to $cos\left(\mathit{A}\right)$, hence one writes $\mathit{W}=\kappa cos\left(\mathit{A}\right)$ [Cannon, 1998, Köhler, 1998].

If knowledge about the value of $\mathit{A}$ is modeled by a Gaussian distribution with mean $\alpha >0$ and standard measurement uncertainty $\mathit{u}\left(\mathit{A}\right)$ (both expressed in radians), then §5.2 (with $\mathit{m}=1$) suggests that knowledge of $\mathit{W}$ should be described approximately by a Gaussian distribution with mean $\kappa cos\left(\alpha \right)$ and standard measurement uncertainty $\kappa \mathit{u}\left(\mathit{A}\right)sin\left(\alpha \right)$.

If, for example, $\alpha =\text{}\pi \text{/3}\phantom{\rule{0.3em}{0ex}}\text{rad}$, $\kappa =\text{2}\phantom{\rule{0.3em}{0ex}}\text{W}$, and $\mathit{u}\left(\mathit{A}\right)=\pi \u2215100\phantom{\rule{0.3em}{0ex}}\text{rad}$, then the approximation to $\mathit{W}$’s distribution that the Delta Method suggests, of a Gaussian distribution with mean $cos\left(\pi \u22153\right)=\text{0.5}\phantom{\rule{0.3em}{0ex}}\text{W}$ and standard deviation $2\left(\pi \u2215100\right)sin\left(\pi \u22153\right)=\text{0.0544}\phantom{\rule{0.3em}{0ex}}\text{W}$, is remarkably accurate (Figure 3).

When the detector is aimed squarely at the target, that is $\alpha =0$, this approximation no longer works because the first derivative of the cosine vanishes at $0$, which is the degenerate case that §5.3 contemplates. In this case, $\mathit{W}$’s standard measurement uncertainty is approximately $\kappa {\mathit{u}}^{2}\left(\mathit{A}\right)\u2215\sqrt{2}$. For $\kappa =\text{2}\phantom{\rule{0.3em}{0ex}}\text{W}$ and $\mathit{u}\left(\mathit{A}\right)=\text{}\pi \text{/100}\phantom{\rule{0.3em}{0ex}}\text{rad}$, this equals $\text{0.0014}\phantom{\rule{0.3em}{0ex}}\text{W}$, which is accurate to the two significant digits shown. However, Figure 3 shows that, in this case, the Delta Method produces a poor approximation.

When $\alpha =0$, the probability density of $\mathit{W}$ is markedly asymmetrical, and the meaning of $\mathit{W}$’s standard deviation is rather different from its meaning when $\alpha >0$. Indeed, when $\alpha =0$ the probability that $\mathit{W}$ should lie within one standard deviation of its expected value is $88\phantom{\rule{0.3em}{0ex}}\%$ approximately.

The Delta Method can be extended to apply to a function of several random variables. Suppose that ${\mathit{X}}_{1}=\left({\mathit{V}}_{1,1}+\cdots +{\mathit{V}}_{{\mathit{m}}_{1},1}\right)$, …, ${\mathit{X}}_{\mathit{n}}=\left({\mathit{V}}_{1,\mathit{n}}+\cdots +{\mathit{V}}_{{\mathit{m}}_{\mathit{n}},\mathit{n}}\right)$ are averages of sets of random variables whose variances are finite. The variables in each set are independent and identically distributed, those in set $1\le \mathit{j}\le \mathit{n}$ having mean ${\mu}_{\mathit{j}}$ and variance ${\sigma}_{\mathit{j}}^{2}$. However, variables in different sets may be dependent, hence the $\left\{{\mathit{X}}_{\mathit{j}}\right\}$ may be dependent, too. Let $\Sigma $ denote the $\mathit{n}\times \mathit{n}$ symmetrical matrix whose element ${\sigma}_{{\mathit{j}}_{1}{\mathit{j}}_{2}}$ $=\mathbb{\mathbb{E}}\left(\right.\left({\mathit{V}}_{\mathit{m},{\mathit{j}}_{1}}-{\mu}_{{\mathit{j}}_{1}}\right)\left({\mathit{V}}_{\mathit{i},{\mathit{j}}_{2}}-{\mu}_{{\mathit{j}}_{2}}\right)\left)\right.$ is the covariance between ${\mathit{X}}_{{\mathit{j}}_{1}}$ and ${\mathit{X}}_{{\mathit{j}}_{2}}$, for $1\le {\mathit{j}}_{1},{\mathit{j}}_{2}\le \mathit{n}$.

Now, consider the random variable $\mathit{Y}=\phi \left({\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}\right)$, where $\phi $ denotes a real-valued function of $\mathit{n}$ variables whose first partial derivatives are continuous and none vanishes at ${\mu}_{1},\dots ,{\mu}_{\mathit{n}}$. If ${\tau}^{2}$ $={\sum}_{{\mathit{j}}_{1}=1}^{\mathit{n}}{\sum}_{{\mathit{j}}_{2}=1}^{\mathit{n}}{\sigma}_{{\mathit{j}}_{1}{\mathit{j}}_{2}}\left(\partial \phi \u2215\partial {\mu}_{{\mathit{j}}_{1}}\right)\left(\mu \right)$ $\left(\partial \phi \u2215\partial {\mu}_{{\mathit{j}}_{2}}\right)\left(\mu \right)$ is finite, then $\sqrt{\mathit{m}}\left(\phi \left(\mathit{Y}\right)-\phi \left({\mu}_{1},\dots ,{\mu}_{\mathit{n}}\right)\right)$ also is approximately Gaussian with mean $0$ and variance ${\tau}^{2}$.

If ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ are uncorrelated and have means ${\mu}_{1},\dots ,{\mu}_{\mathit{n}}$ and standard deviations ${\sigma}_{1},\dots ,{\sigma}_{\mathit{n}}$, then the Delta Method approximation reduces to a well-known formula first presented by Gauss [Gauss, 1823, §18, Problema]: $\mathbb{\mathbb{V}}\left(\mathit{Y}\right)\approx {\mathit{c}}_{1}^{2}{\sigma}_{1}^{2}+\cdots +{\mathit{c}}_{\mathit{n}}^{2}{\sigma}_{\mathit{n}}^{2}$, where the sensitivity coefficient ${\mathit{c}}_{\mathit{j}}=\partial \phi \left({\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)\u2215\partial {\mathit{x}}_{\mathit{j}}$ is the value at $\left({\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)$ of the $\mathit{j}$th partial derivative of $\phi $ with respect to ${\mathit{x}}_{\mathit{j}}$.

If a beam of monochromatic light of power
${\mathit{I}}_{0}$ (W) travels a
path of length $\mathit{L}$
(m) through a solution containing a solute whose molar absorptivity for that light
is $\mathit{E}$
(L mol^{-1} m^{-1}), and whose molar concentration is
$\mathit{C}$
(mol L^{-1}), then the beam’s power is reduced to
$\mathit{I}$ (W) such that
$\mathit{I}={\mathit{I}}_{0}{10}^{\mathit{E}\mathit{L}\mathit{C}}$. Application of Gauss’s
formula (§5.5) to $\mathit{C}={log}_{10}\left({\mathit{I}}_{0}\u2215\mathit{I}\right)\u2215\left(\mathit{E}\mathit{L}\right)$
yields:

$$\mathbb{\mathbb{V}}\left(\mathit{C}\right)\approx \frac{\frac{{\sigma}_{{\mathit{I}}_{0}}^{2}}{{\mathit{I}}_{0}^{2}}+\frac{{\sigma}_{\mathit{I}}^{2}}{{\mathit{I}}^{2}}}{{\left(\mathit{E}\mathit{L}log10\right)}^{2}}+\left(\right.\frac{{\sigma}_{\mathit{E}}^{2}}{{\mathit{E}}^{2}}+\frac{{\sigma}_{\mathit{L}}^{2}}{{\mathit{L}}^{2}}\left)\right.{log}_{10}^{2}\left(\mathit{I}\u2215{\mathit{I}}_{0}\right)$$ |

Darcy’s law relates the dynamic viscosity $\mathit{H}$ of a fluid to the volumetric rate of discharge $\mathit{Q}$ (volume per unit of time) when the fluid flows through a permeable cylindrical medium of cross-section $\mathit{A}$ and intrinsic permeability $\mathit{K}$ under a pressure drop of $\Delta $ along a length $\mathit{L}$, as follows: $\mathit{H}=\mathit{K}\mathit{A}\Delta \u2215\left(\mathit{Q}\mathit{L}\right)$.

To compute an approximation to the standard deviation of $\mathit{H}$, one may use the formula from §5.5 directly, or first take the logarithm of both sides, which linearizes the relationship, $log\left(\mathit{H}\right)=log\left(\mathit{K}\right)+log\left(\mathit{A}\right)+log\left(\Delta \right)-log\left(\mathit{Q}\right)-log\left(\mathit{L}\right)$. Applied to these logarithms, the formula from §5.5 is exact. The approximation is then done for each term separately, using the univariate Delta Method.

Since $\mathbb{\mathbb{V}}\left(log\left(\eta \right)\right)\approx \mathbb{\mathbb{V}}\left(\eta \right)\u2215{\eta}^{2}$, and similarly for the other logarithmic terms, $\mathbb{\mathbb{V}}\left(\eta \right)\u2215{\eta}^{2}\approx \mathbb{\mathbb{V}}\left(\kappa \right)\u2215{\kappa}^{2}\cdots +\mathbb{\mathbb{V}}\left(\mathit{L}\right)\u2215{\mathit{L}}^{2}$. In other words, the square of the coefficient of variation of $\eta $ is approximately equal to the sum of the squares of the variation coefficients of the other variables.

The Monte Carlo method offers several important advantages over the Delta Method described in §5.2 and §5.5: (i) it can produce as many correct significant digits in its results as may be required; (ii) it does not involve the computation of derivatives, either analytically or numerically; (iii) it is applicable in many situations where the Delta Method is not; (iv) it provides a picture of the whole probability distribution of a function of several random variables, not just an approximation to it, or to its mean and standard deviation.

The Monte Carlo method in general dates back to the middle of the twentieth century [Metropolis and Ulam, 1949, Metropolis et al., 1953]. A variant used in mathematical statistics is known as the parametric bootstrap [Efron and Tibshirani, 1993]. This involves using random draws from a (possibly multivariate) probability distribution whose parameters have been replaced by estimates thereof (for example, means of posterior probability distributions, §6) to ascertain the probability distribution of a function of one (or more) random variables. Morgan and Henrion [1992] and Joint Committee for Guides in Metrology [2008b] describe how it may be employed to evaluate measurement uncertainty, and provide illustrative examples.

The procedure comprises the following steps:

- MC1
- Define the joint probability distribution of ${\mathit{X}}_{1}$, …, ${\mathit{X}}_{\mathit{n}}$.
- MC2
- Choose a suitably large positive integer $\mathit{K}$ and draw a sample of size $\mathit{K}$ from this joint distribution to obtain $\left({\mathit{x}}_{11},\dots ,{\mathit{x}}_{\mathit{n}1}\right)$, …, $\left({\mathit{x}}_{1\mathit{K}},\dots ,{\mathit{x}}_{\mathit{n}\mathit{K}}\right)$. (If ${\mathit{X}}_{1}$, …, ${\mathit{X}}_{\mathit{n}}$ happen to be independent, then this amounts to drawing a sample of size $\mathit{K}$ from the distribution of each of them separately.)
- MC3
- Compute ${\mathit{y}}_{1}=\phi \left({\mathit{x}}_{11},\dots ,{\mathit{x}}_{\mathit{n}1}\right)$, …, ${\mathit{y}}_{\mathit{K}}=\phi \left({\mathit{x}}_{1\mathit{K}},\dots ,{\mathit{x}}_{\mathit{n}\mathit{K}}\right)$, which are a sample from $\mathit{Y}$’s distribution.
- MC4
- Summarize this sample in one or more of these different ways:
- MC4.a — Probability Density
- The most inclusive summarization is in the form of an estimate of $\mathit{Y}$’s probability density function: this may be either a simple histogram, or a kernel density estimate [Silverman, 1986].
- MC4.b — Mean and Standard Deviation
- The mean and standard deviation of $\mathit{Y}$ are estimated by the mean and the standard deviation of $\left\{{\mathit{y}}_{1},\dots ,{\mathit{y}}_{\mathit{K}}\right\}$. (To ascertain the number of significant digits in this mean and standard deviation, hence to decide whether $\mathit{K}$ is large enough for the intended purpose, or should be increased, one may employ either the adaptive procedure explained in the Supplement 1 to the GUM [Joint Committee for Guides in Metrology, 2008b, 7.9], or resort to the non-parametric statistical bootstrap or to other resampling methods [Davison and Hinkley, 1997].)
- MC4.c — Probability Interval
- If ${\mathit{y}}_{\left(1\right)}\le {\mathit{y}}_{\left(2\right)}\le \cdots \le {\mathit{y}}_{\left(\mathit{K}\right)}$ denote the result of ordering ${\mathit{y}}_{1},\dots ,{\mathit{y}}_{\mathit{K}}$ from smallest to largest, then the interval $\left({\mathit{y}}_{\left(\mathit{K}\alpha \u22152\right)},{\mathit{y}}_{\left(\mathit{K}\left(1-\alpha \u22152\right)\right)}\right)$ includes $\mathit{Y}$’s true value with probability $1-\alpha $. (Since $\mathit{K}\alpha \u22152$ and $\mathit{K}\left(1-\alpha \u22152\right)$ need not be integers, the end-points of this coverage interval may be calculated by interpolation of adjacent ${\mathit{y}}_{\left(\mathit{i}\right)}$s.)

The radius $\mathit{R}$
and the height $\mathit{H}$
of a cylinder are values of independent random variables with exponential probability
distributions with mean 1 m. To characterize the probability distribution of its volume
$\mathit{V}=\pi {\mathit{R}}^{2}\mathit{H}$, draw a sample of
size ${10}^{7}$ from the joint
distribution of $\mathit{R}$
and $\mathit{H}$,
and compute the volume corresponding to each pair of sampled values
$\left(\mathit{r},\mathit{h}\right)$
to obtain a sample of the same size from the distribution of
$\mathit{V}$. The
average and standard deviation of these values are 6.3 m^{3} and 21 m^{3}, and they are
estimates (whose two most significant digits are exact) of the mean and standard
deviation of $\mathit{V}$.
Figure 4 depicts an estimate of the corresponding probability density.

Suppose that $\mathit{X}$ is a random variable with a continuous distribution and values in $\mathcal{\mathcal{X}}$, with probability distribution function ${\mathit{P}}_{\mathit{X}}$ and probability density function ${\mathit{p}}_{\mathit{X}}$, and consider the random variable $\mathit{Y}=\phi \left(\mathit{X}\right)$ where $\phi $ denotes a real-valued function of a real variable. Let $\mathcal{\mathcal{Y}}$ denote the set where $\mathit{Y}$ takes its values, and let ${\mathit{P}}_{\mathit{Y}}$ and ${\mathit{p}}_{\mathit{Y}}$ denote $\mathit{Y}$’s probability distribution function and probability density function, respectively. In these circumstances [Casella and Berger, 2002, Chapter 2].

- If $\phi $ is increasing on $\mathcal{\mathcal{X}}$ and $\psi $ denotes its inverse, then ${\mathit{P}}_{\mathit{Y}}\left(\mathit{y}\right)=Pr\left(\mathit{Y}\le \mathit{y}\right)=Pr\left(\mathit{X}\le \psi \left(\mathit{y}\right)\right)={\mathit{P}}_{\mathit{X}}\left[\psi \left(\mathit{y}\right)\right]$ for $\mathit{y}\in \mathcal{\mathcal{Y}}$; and if $\phi $ is decreasing, then ${\mathit{P}}_{\mathit{Y}}\left(\mathit{y}\right)=1-{\mathit{P}}_{\mathit{X}}\left[\psi \left(\mathit{y}\right)\right]$.
- If $\phi $ is either increasing or decreasing on $\mathcal{\mathcal{X}}$ (but not both), and its inverse $\psi $ has a continuous first derivative $\stackrel{\u0307}{\psi}$, then ${\mathit{p}}_{\mathit{Y}}\left(\mathit{y}\right)={\mathit{p}}_{\mathit{X}}\left[\psi \left(\mathit{y}\right)\right]\left|\stackrel{\u0307}{\psi}\left(\mathit{y}\right)\right|$, for $\mathit{y}\in \mathcal{\mathcal{Y}}$, where $\left|\stackrel{\u0307}{\psi}\left(\mathit{y}\right)\right|$ denotes the absolute value of the derivative of $\psi $ at $\mathit{y}$.

A horizontal beam of light emerges from a tiny hole in a wall and travels along a 1 m long path at right angles to the wall, towards a flat mirror that oscillates freely around a vertical axis. When the mirror’s surface normal makes an angle $\mathit{A}$ with the beam, its reflection hits the wall at distance $\mathit{D}=tan\left(\mathit{A}\right)$ from the hole (positive to the right of the hole and negative to the left). If $\mathit{A}$ is uniformly (or, rectangularly) distributed between $\text{-}\pi \text{/2}\phantom{\rule{0.3em}{0ex}}\text{rad}$ and $\text{}\pi \text{/2}\phantom{\rule{0.3em}{0ex}}\text{rad}$, then ${\mathit{P}}_{\mathit{D}}\left(\mathit{d}\right)=Pr\left(\mathit{D}\le \mathit{d}\right)=Pr\left(\mathit{A}\le arctan\left(\mathit{d}\right)\right)=\left(arctan\left(\mathit{d}\right)+\pi \u22152\right)\u2215\pi $, and $\mathit{D}$’s probability density is ${\mathit{p}}_{\mathit{D}}$ such that ${\mathit{p}}_{\mathit{D}}\left(\mathit{d}\right)=1\u2215\left[\pi \left(1+{\mathit{d}}^{2}\right)\right]$ for $-\infty <\mathit{d}<\infty $. As it turns out, both the mean and the standard deviation of $\mathit{D}$ are infinite [Feller, 1971, Page 51].

Suppose that ${\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}$ are random variables and consider ${\mathit{Y}}_{\mathit{j}}={\phi}_{\mathit{j}}\left({\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}\right)$ for $\mathit{j}=1,\dots ,\mathit{n}$, and where ${\phi}_{1},\dots ,{\phi}_{\mathit{n}}$ are real-valued functions of $\mathit{n}$ real variables each.

Suppose also that (i) the vector $\mathit{X}=\left({\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}\right)$ takes values in an open subset $\mathcal{\mathcal{X}}$ of $\mathit{n}$-dimensional Euclidean space, and has a continuous joint probability distribution with probability density function ${\mathit{p}}_{\mathit{X}}$; (ii) the vector-valued function $\phi =\left({\phi}_{1},\dots ,{\phi}_{\mathit{n}}\right)$ is invertible, and the inverse $\psi =\left({\psi}_{1},\dots ,{\psi}_{\mathit{n}}\right)$ has a Jacobian determinant ${\mathit{J}}_{\psi}$ that does not vanish on $\mathcal{\mathcal{Y}}$.

- CV1
- Solve the $\mathit{n}$ equations ${\mathit{y}}_{1}={\phi}_{1}\left({\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)$, …, ${\mathit{y}}_{\mathit{n}}={\phi}_{\mathit{n}}\left({\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)$, for ${\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}$, to obtain the inverse transformation such that ${\mathit{x}}_{1}={\psi}_{1}\left({\mathit{y}}_{1},\dots ,{\mathit{y}}_{\mathit{n}}\right)$, …, ${\mathit{x}}_{\mathit{n}}={\psi}_{\mathit{n}}\left({\mathit{y}}_{1},\dots ,{\mathit{y}}_{\mathit{n}}\right)$.
- CV2
- Find ${\stackrel{\u0307}{\psi}}_{\mathit{i}\mathit{j}}$, the partial
derivative of ${\psi}_{\mathit{i}}$ with
respect to its $\mathit{j}$th
argument, for $\mathit{i},\mathit{j}=1,\dots ,\mathit{n}$,
and compute the Jacobian determinant of the inverse transformation at
$\mathit{y}$
$=\left({\mathit{y}}_{1},\dots ,{\mathit{y}}_{\mathit{n}}\right)$:
$${\mathit{J}}_{\psi}\left(\mathit{y}\right)=det\left[\begin{array}{cccc}\hfill {\stackrel{\u0307}{\psi}}_{11}\left(\mathit{y}\right)\hfill & \hfill {\stackrel{\u0307}{\psi}}_{12}\left(\mathit{y}\right)\hfill & \hfill \dots \hfill & \hfill {\stackrel{\u0307}{\psi}}_{1\mathit{n}}\left(\mathit{y}\right)\hfill \\ \hfill {\stackrel{\u0307}{\psi}}_{21}\left(\mathit{y}\right)\hfill & \hfill {\stackrel{\u0307}{\psi}}_{22}\left(\mathit{y}\right)\hfill & \hfill \dots \hfill & \hfill {\stackrel{\u0307}{\psi}}_{2\mathit{n}}\left(\mathit{y}\right)\hfill \\ \hfill \vdots \hfill & \hfill \vdots \hfill & \hfill \ddots \hfill & \hfill \vdots \hfill \\ \hfill {\stackrel{\u0307}{\psi}}_{\mathit{n}1}\left(\mathit{y}\right)\hfill & \hfill {\stackrel{\u0307}{\psi}}_{\mathit{n}2}\left(\mathit{y}\right)\hfill & \hfill \dots \hfill & \hfill {\stackrel{\u0307}{\psi}}_{\mathit{n}\mathit{n}}\left(\mathit{y}\right)\hfill \\ \hfill \hfill \end{array}\right]$$ - CV3
- The density the joint probability distribution of the random vector
$\mathit{Y}$ is
${\mathit{p}}_{\mathit{Y}}$
such that
$${\mathit{p}}_{\mathit{Y}}\left(\mathit{y}\right)={\mathit{p}}_{\mathit{X}}\left[\psi \left(\mathit{y}\right)\right]\left|{\mathit{J}}_{\psi}\left(\mathit{y}\right)\right|.$$ (2) Note that ${\mathit{J}}_{\psi}\left(\mathit{y}\right)$ is a scalar, and $\left|{\mathit{J}}_{\psi}\left(\mathit{y}\right)\right|$ denotes its absolute value.

- CV4
- The probability density of ${\mathit{Y}}_{1}$ is ${\mathit{p}}_{{\mathit{Y}}_{1}}\left({\mathit{y}}_{1}\right)=\int \dots \int \mathit{g}\left({\mathit{y}}_{1},{\mathit{y}}_{2},\dots ,{\mathit{y}}_{\mathit{n}}\right)\text{d}{\mathit{y}}_{2}\dots \text{d}{\mathit{y}}_{\mathit{n}}$, where the $\mathit{n}-1$ integrals are over the ranges of ${\mathit{Y}}_{2},\dots ,{\mathit{Y}}_{\mathit{n}}$.

Suppose that $\mathit{U}$ and $\mathit{V}$ are independent, Gaussian random variables with mean 0 and variance 1, and let $\mathit{S}=\mathit{a}\mathit{U}+\mathit{b}\mathit{V}$, and $\mathit{T}=\mathit{b}\mathit{U}-\mathit{a}\mathit{V}$, for given real numbers $\mathit{a}$ and $\mathit{b}$. The inverse transformation maps $\left(\mathit{s},\mathit{t}\right)$ onto $\left(\right.\left(\mathit{a}\mathit{s}+\mathit{b}\mathit{t}\right)\u2215\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right),\left(\mathit{b}\mathit{s}-\mathit{a}\mathit{t}\right)\u2215\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)\left)\right.$, and has Jacobian determinant $det\left[\begin{array}{cc}\mathit{a}& \mathit{b}\\ \mathit{b}& -\mathit{a}\end{array}\right]\u2215\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)$ whose absolute value is $1\u2215\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)$. Since the density of the joint probability distribution of $\mathit{U}$ and $\mathit{V}$ is ${\mathit{p}}_{\mathit{U},\mathit{V}}\left(\mathit{u},\mathit{v}\right)=exp\left(-\left({\mathit{u}}^{2}+{\mathit{v}}^{2}\right)\u22152\right)\u2215\left(2\pi \right)$, application of the multivariate change-of-variable formula yields

$${\mathit{p}}_{\mathit{S},\mathit{T}}\left(\mathit{s},\mathit{t}\right)=\frac{exp\left\{-\frac{{\mathit{s}}^{2}}{2\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)}\right\}}{\sqrt{2\pi \left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)}}\frac{exp\left\{-\frac{{\mathit{t}}^{2}}{2\left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)}\right\}}{\sqrt{2\pi \left({\mathit{a}}^{2}+{\mathit{b}}^{2}\right)}}.$$ |

But this means that $\mathit{S}$ and $\mathit{T}$ also are independent and Gaussian with mean 0 and variance ${\mathit{a}}^{2}+{\mathit{b}}^{2}$. On the one hand, this result is surprising because $\mathit{S}$ and $\mathit{T}$ both are functions of the same random variables. On the other hand it is hardly surprising because the transformation amounts to a rotation of the coordinate axes, followed by a global dilation. Since the joint distribution of $\mathit{U}$ and $\mathit{V}$ is circularly symmetric relative to $\left(0,0\right)$, so will the joint distribution of $\mathit{S}$ and $\mathit{T}$ be, which implies independence and the same functional form for the density, up to a difference in scale.

To compute the probability density of the ratio $\mathit{R}=\mathit{X}\u2215\mathit{Y}$ of two independent and exponentially distributed random variables $\mathit{X}$ and $\mathit{Y}$ with mean $1\u2215\lambda $, define the function $\phi $ such that $\phi \left(\mathit{x},\mathit{y}\right)=\left(\mathit{x}\u2215\mathit{y},\mathit{y}\right)$, whose inverse is $\psi \left(\mathit{r},\mathit{s}\right)=\left(\mathit{r}\mathit{s},\mathit{s}\right)$, with Jacobian determinant ${\mathit{J}}_{\psi}\left(\mathit{r},\mathit{s}\right)=det\left[\begin{array}{cc}\mathit{s}& \mathit{r}\\ 0& 1\end{array}\right]=\mathit{s}>0$. The multivariate change-of-variable formula then yields ${\mathit{p}}_{\mathit{R},\mathit{S}}\left(\mathit{r},\mathit{s}\right)=\mathit{s}{\lambda}^{2}exp\left[\right.-\lambda \left(1+\mathit{r}\right)\mathit{s}\left]\right.$ for the density of the joint distribution of $\mathit{R}=\mathit{X}\u2215\mathit{Y}$ and $\mathit{S}=\mathit{Y}$. The (marginal) density of $\mathit{R}$ is ${\mathit{p}}_{\mathit{R}}\left(\mathit{r}\right)={\int}_{0}^{\infty}{\mathit{p}}_{\mathit{R},\mathit{S}}\left(\mathit{r},\mathit{s}\right)\text{d}\mathit{s}=1\u2215{\left(1+\mathit{r}\right)}^{2}$, for $\mathit{r}>0$, being 0 otherwise.

${\mathit{p}}_{\mathit{R}}$ indeed is a probability density function because it is a non-negative function and ${\int}_{0}^{\infty}\text{d}\mathit{r}\u2215{\left(1+\mathit{r}\right)}^{2}=1-\underset{\mathit{r}\to \infty}{lim}1\u2215\left(1+\mathit{r}\right)=1$. However, since ${\int}_{0}^{\infty}\mathit{r}\text{d}\mathit{r}\u2215{\left(1+\mathit{r}\right)}^{2}=\underset{\mathit{r}\to \infty}{lim}log\left(1+\mathit{r}\right)=\infty $, neither the mean nor the variance of $\mathit{R}$ is finite. The Delta Method, however, would have suggested that $\mathbb{\mathbb{E}}\left(\mathit{X}\u2215\mathit{Y}\right)\approx 1$ and that the coefficient of variation (ratio of the standard deviation to the mean) of $\mathit{X}\u2215\mathit{Y}$ is $\sqrt{2}$ approximately.

The inverse of the transformation that maps the polar coordinates $\left(\mathit{r},\alpha \right)$ of a point in the Euclidean plane to its Cartesian coordinates $\left(\mathit{x},\mathit{y}\right)$ is $\psi $ such that $\psi \left(\mathit{r},\mathit{a}\right)=\left(\mathit{r}cos\mathit{a},\mathit{r}sin\mathit{a}\right)$ for $\mathit{r}>0$ and $0<\mathit{a}<2\pi $, with Jacobian determinant ${\mathit{J}}_{\psi}\left(\mathit{r},\alpha \right)=\mathit{r}$. If $\mathit{X}$ and $\mathit{Y}$ are independent Gaussian random variables with mean 0 and variance 1, then the probability density of $\left(\mathit{R},\mathit{A}\right)$ is ${\mathit{p}}_{\mathit{R},\mathit{A}}\left(\mathit{r},\alpha \right)=\mathit{r}exp\left(-{\mathit{r}}^{2}\u22152\right)\u2215\left(2\pi \right)$. Since ${\int}_{0}^{\infty}\mathit{r}exp\left(-{\mathit{r}}^{2}\u22152\right)=1$, it follows that $\mathit{R}$ and $\mathit{A}$ are independent, the former having a Rayleigh distribution with mean $\sqrt{\pi \u22152}$, the latter a uniform distribution between $0$ and $2\pi $.

The statistical inferences we are primarily interested in are probabilistic statements about the unknown value of a quantity, produced by application of a statistical method. In the example of §6.6, one of the inferences is this statement: the probability is 95% that the difference between the numbers of hours gained with the two soporifics is between 0.7 h and 2.5 h.

Another common inference is an estimate of the value of a quantity, which must be qualified with an assessment of the associated uncertainty. In the example of §6.8, a typical inference of this kind would be this: the difference in mean levels of thyroxine in the serum of two groups of children diagnosed with hypothyroidism is estimated as $\text{14}\phantom{\rule{0.3em}{0ex}}\text{nmol/L}$ with standard uncertainty $\text{18}\phantom{\rule{0.3em}{0ex}}\text{nmol/L}$ (Figure 6).

In our treatment of this example, the inference is based entirely on a small set of empirical data, and on a particular choice of statistical model used to describe the dispersion of the data, and to characterize the fact that no knowledge other than the data was brought into play.

Statistical methods different from the one we used could have been employed: some of these would produce the same result (in particular those illustrated when this dataset was first described [Student, 1908, Fisher, 1973]), while others would have produced different results.

Even when the result is the same, it may be variously interpreted:

- For some that statement means that if the same sampling and study method is used repeatedly, and each time the resulting dataset is modeled and analyzed in the same way to produce an interval like the one above, then about 95 % of the resulting intervals will include the true difference sought — with no guarantee or implication that the interval that was obtained is one of these;
- For others (among whom we stand) that statement expresses the degree of belief one is entitled to have about the true difference lying between 0.7 h and 2.5 h specifically, in light of all the relevant information in hand.

Bayesian inference [Bernardo and Smith, 2000, Lindley, 2006, Robert, 2007] is a class of statistical procedures that serve to blend preexisting information about the value of a quantity with fresh information in empirical data.

The defining traits of a Bayesian procedure are these:

- (i)
- All quantity values that are the objects of interest but are accessible to direct observation (non-observables) are modeled as values of non-observable random variables whose (prior, or a priori) distributions encode and convey states of incomplete knowledge about those values;
- (ii)
- The empirical data (observables) are modeled as realized values of random variables whose probability distributions depend on those objects of interest;
- (iii)
- Preexisting information about those objects of interest is updated in light of the fresh empirical data by application of Bayes rule, and the results are encapsulated in a (posterior, or a posteriori) probability distribution;
- (iv)
- Selected aspects of this distribution are then abstracted from it and used to characterize the objects of interest and to describe the state of knowledge about them.

Let $\theta $ denote the value of the quantity of interest, which we model as realized value of a random variable $\Theta $ with probability density function ${\mathit{p}}_{\Theta}$ that encodes the state of knowledge about $\theta $ prior to obtaining fresh data, and which must be defined even if there is no prior knowledge.

Defining such ${\mathit{p}}_{\Theta}$ often is a challenging task. If in fact there exists substantial prior knowledge about $\theta $, then it needs to be elicited from experts in the matter and encapsulated in the form of a particular probability density: Garthwaite et al. [2005] review how this may be done. For example, when measuring the mass fraction of titanium in a mineral specimen, then knowledge of the species (ilmenite, titanite, rutile, etc.) of the specimen is highly informative about that mass fraction. Familiarity with the process of analytical chemistry employed to make the measurement may indicate the dispersion of values to be expected.

In some cases, essentially no prior knowledge exists about $\theta $, or none is deemed reliable enough to be taken into account. In such cases, a so-called non-informative prior distribution needs to be produced and assigned to $\Theta $ that reflects this state of affairs: if $\theta $ is univariate (that is, a single number), then the rules developed by Jeffreys [1961] often prove satisfactory; if $\theta $ is multivariate (that is, a numerical vector), then the so-called reference prior distributions are recommended [Bernardo and Smith, 2007] (these reduce to Jeffreys’s in the univariate case).

These rules often produce a ${\mathit{p}}_{\Theta}$ that is improper, in the sense that ${\int}_{\mathcal{\mathscr{H}}}\mathit{p}\left(\theta \right)\text{d}\theta $ diverges to infinity, where $\mathcal{\mathscr{H}}$ denotes the range of $\Theta $. (If $\Theta $ should have a discrete distribution then this integral is replaced by a sum.) Fortunately, once used in Bayes Rule (§3.5 and §6.4), improper priors often lead to proper posterior probability distributions.

The empirical data $\mathit{x}$ (which may be a single number, a numerical vector, or a data structure of still greater complexity) are modeled as realized values of a random variable $\mathit{X}$ whose probability density describes the corresponding dispersion of values.

This density must depend on $\theta $, which is another way of saying that the data are informative about $\theta $ (otherwise there would be nothing to be gained by observing them). In fact, this is the density of the conditional probability distribution of $\mathit{X}$ given that $\Theta =\theta $. Choosing a specific functional form for it generally is a non-trivial exercise: it involves defining a statistical model that correctly captures the dispersion of values likely to be obtained in the experiment that produces them.

Once the data $\mathit{x}$ are in hand, ${\mathit{p}}_{\mathit{X}|\Theta}\left(\mathit{x}|\theta \right)$ becomes a function of $\theta $ alone, being largest for values of $\theta $ that make the data appear most likely. As such, it still is non-negative, but its integral (or sum, if $\mathit{X}$’s distribution should be discrete) over the range of $\Theta $, need not be 1.

Suppose that both $\mathit{X}$ given that $\Theta =\theta $ and $\Theta $ have continuous distributions with densities ${\mathit{p}}_{\mathit{X}|\Theta}$ and ${\mathit{p}}_{\Theta}$. In these circumstances, Bayes rule becomes

The function ${\mathit{p}}_{\Theta |\mathit{X}}$, which is defined over the range of $\Theta $ for each fixed value $\mathit{x}$, is the density of the posterior distribution of the value of the quantity of interest $\Theta $ given the data.

In some cases this can be computed in closed form, in many others it cannot. In all cases it is possible to obtain a sample from this posterior distribution by application of a procedure known as Markov Chain Monte Carlo (MCMC) [Gelman et al., 2003]. This sample can then be summarized as described in §5.8.

Once infected by influenza A virus, an epithelial cell of the upper respiratory tract releases $\theta $ virions on average, which may then go on to infect other cells. This number $\theta $ depends on the volume of the cell, and we will treat it as realized value of a non-observable random variable with an exponential distribution whose expected value, $1\u2215\gamma $ for some $0<\gamma <1$, is known. Given $\theta $, the actual number of virions that are released is $\mathit{x}$, and this is like a realized value of a Poisson random variable with mean $\theta $.

Suppose that the prior density is ${\mathit{p}}_{\Theta}\left(\theta \right)=\gamma exp\left(-\gamma \theta \right)$, and the likelihood function is ${\mathit{L}}_{\mathit{x}}\left(\theta \right)={\mathit{p}}_{\mathit{X}|\Theta}\left(\mathit{x}|\theta \right)={\theta}^{\mathit{x}}exp\left(-\theta \right)\u2215\mathit{x}!$, for $\theta >0$. The posterior distribution of $\Theta $ given $\mathit{x}$ belongs to the gamma family, and has expected value $\left(\mathit{x}+1\right)\u2215\left(\gamma +1\right)$, variance $\left(\mathit{x}+1\right)\u2215{\left(\gamma +1\right)}^{2}$, and density

$${\mathit{p}}_{\Theta |\mathit{X}}\left(\theta |\mathit{x}\right)=\frac{\frac{{\theta}^{\mathit{x}}{\mathit{e}}^{-\theta}}{\mathit{x}!}\gamma {\mathit{e}}^{-\gamma \theta}}{{\int}_{0}^{\infty}\frac{{\mathit{s}}^{\mathit{x}}{\mathit{e}}^{-\mathit{s}}}{\mathit{x}!}\gamma {\mathit{e}}^{-\gamma \mathit{s}}\text{d}\mathit{s}}=\frac{{\left(\gamma +1\right)}^{\mathit{x}+1}{\theta}^{\mathit{x}}{\mathit{e}}^{-\theta \left(\gamma +1\right)}}{\mathit{x}!}.$$ |

The differences between the numbers of additional hours of sleep that ten patients gained when using two soporific drugs, described in examples given by Student [1908] and [Fisher, 1973, §24], were 1.2 h, 2.4 h, 1.3 h, 1.3 h, 0.0 h, 1.0 h, 1.8 h, 0.8 h, 4.6 h, and 1.4 h.

Suppose that, given $\mu $ and $\sigma $, these are realized values of independent Gaussian random variables with mean $\mu $ and variance ${\sigma}^{2}$. Let $\overline{\mathit{x}}=\text{1.58}\phantom{\rule{0.3em}{0ex}}\text{h}$ denote their average, and ${\mathit{s}}^{2}=1.51\phantom{\rule{0.3em}{0ex}}{\text{h}}^{2}$ denote the sum of their squared deviations from $\overline{\mathit{x}}$ divided by 9. In these circumstances, the likelihood function is ${\mathit{L}}_{\overline{\mathit{x}},{\mathit{s}}^{2}}\left(\mu ,{\sigma}^{2}\right)={\left(2\pi {\sigma}^{2}\right)}^{-\mathit{n}\u22152}exp\left\{-\left[\mathit{n}{\left(\overline{\mathit{x}}-\mu \right)}^{2}+\left(\mathit{n}-1\right){\mathit{s}}^{2}\right]\u2215\left(2{\sigma}^{2}\right)\right\}$.

Assume, in addition, that $\mu $ and $\sigma $ are realized values of non-observable random variables $\mathit{M}$ and $\Sigma $ that are independent a priori and such that $\mathit{M}$ and $log\Sigma $ are uniformly distributed between $-\infty $ and $+\infty $ (both improper prior distributions). Then, given $\overline{\mathit{x}}$ and $\mathit{s}$, $\left(\mu -\overline{\mathit{x}}\right)\u2215\left(\mathit{s}\u2215\sqrt{\mathit{n}}\right)$ is like a realized value of a random variable with a Student’s $\mathit{t}$ distribution with $\mathit{n}-1=9$ degrees of freedom, and $\left(\mathit{n}-1\right){\mathit{s}}^{2}\u2215{\sigma}^{2}$ is like a realized value of a random variable with a chi-squared distribution with $\mathit{n}-1$ degrees of freedom [Box and Tiao, 1973, Theorem 2.2.1].

Therefore, the expected value of the posterior distribution of the mean difference of hours of sleep gained is $\tilde{\mu}=\text{1.58}\phantom{\rule{0.3em}{0ex}}\text{h}$, and the standard deviation is $\left(\right.\mathit{s}\u2215\sqrt{\mathit{n}}\left)\right.\left(\right.\sqrt{\left(\mathit{n}-1\right)\u2215\left(\mathit{n}-3\right)}\left)\right.=\text{0.441}\phantom{\rule{0.3em}{0ex}}\text{h}$. A 95 % probability interval for $\mu $ ranges from 0.7 h to 2.5 h, and a similar one for $\sigma $ ranges from 0.8 h to 2.2 h.

Suppose that ${\mathit{M}}_{\overline{\mathit{x}},\mathit{s}}$ and ${\Sigma}_{\overline{\mathit{x}},\mathit{s}}$ are the counterparts of $\mathit{M}$ and $\Sigma $ once the information in the data has been taken into account: that is, their probability distribution is the joint (or, bivariate) posterior probability distribution given the data. Even though $\mathit{M}$ and $\Sigma $ were assumed to be independent a priori, ${\mathit{M}}_{\overline{\mathit{x}},\mathit{s}}$ and ${\Sigma}_{\overline{\mathit{x}},\mathit{s}}$ turn out to be dependent a posteriori (that is, given the data), but their correlation is zero [Lindley, 1965b, §5.4].

A major hurricane has category 3, 4, or 5 on the Saffir-Simpson Hurricane Scale
[Simpson, 1974]: its central pressure is no more than 945 mbar (94 500 Pa), it has winds
of at least $111\phantom{\rule{0.3em}{0ex}}\text{mile}\u2215\text{hour}$
(49.6 m s^{-1}), generates sea surges of 9 feet (2.7 m) or greater, and has the
potential to cause extensive damage.

The numbers of major hurricanes that struck the U.S. mainland directly, in each decade starting with 1851–1860 and ending with 1991–2010, are: 6, 1, 7, 5, 8, 4, 7, 5, 8, 10, 8, 6, 4, 4, 5, 7 [Blake et al., 2011]. Let $\mathit{n}=16$ denote the number of decades, ${\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}$ denote the corresponding counts, and $\mathit{s}={\mathit{x}}_{1}+\cdots +{\mathit{x}}_{\mathit{n}}$. Suppose that one wishes to predict $\mathit{y}$, the number of such hurricanes in the decade 2011–2020.

Assume that the mean number of such hurricanes per decade will have remained constant between 1851 and 2010 (certainly a questionable assumption), with unknown value $\lambda $, and that, conditionally on this value, ${\mathit{x}}_{1}$, …, ${\mathit{x}}_{\mathit{n}}$, and $\mathit{y}$ are realized values of independent Poisson random variables ${\mathit{X}}_{1}$, …, ${\mathit{X}}_{\mathit{n}}$ (observable), $\mathit{Y}$ (non-observable), all with mean value $\lambda $: their common probability density is ${\mathit{p}}_{\mathit{X}|\Lambda}\left(\mathit{k}|\lambda \right)={\lambda}^{\mathit{k}}exp\left(-\lambda \right)\u2215\mathit{k}!$ for $\mathit{k}=0,1,2,\dots \phantom{\rule{0.3em}{0ex}}$. This model is commonly used for phenomena that result from the cumulative effect of many improbable events [Feller, 1968, XI.6b].

Even though the goal is to predict $\mathit{Y}$, the fact that there is no a priori knowledge about $\lambda $ other than that it must be positive, requires that this be modeled as the (non-observable) value of a random variable $\Lambda $ whose probability distribution must reflect this ignorance. (According to the Bayesian paradigm, all states of knowledge, even complete ignorance, have to be modeled using probability distributions.)

If the prior distribution chosen for $\Lambda $ is the reference prior distribution [Berger, 2006, Bernardo, 1979, Bernardo and Smith, 2000], then the value of its probability density ${\mathit{p}}_{\Lambda}$ at $\lambda $ should be proportional to $1\u2215\sqrt{\lambda}$ [Bernardo and Smith, 2000, A.2], an improper prior probability density. However, the corresponding posterior distribution for $\Lambda $ is proper, in fact it is a gamma distribution with expected value $\left(\mathit{s}+\text{\xbd}\right)\u2215\mathit{n}$ and probability density function ${\mathit{p}}_{\Lambda |{\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}}$ such that

However, what is needed for the aforementioned prediction is the conditional distribution of $\mathit{Y}$ given the observed counts: the so-called predictive distribution [Schervish, 1995, Page 18]. If $\pi $ denotes the corresponding density, then

$$\begin{array}{llll}\hfill \pi \left(\mathit{y}\phantom{\rule{0.3em}{0ex}}|\phantom{\rule{0.3em}{0ex}}{\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)& ={\int}_{0}^{+\infty}\mathit{g}\left(\mathit{y}\phantom{\rule{0.3em}{0ex}}|\phantom{\rule{0.3em}{0ex}}\lambda \right){\mathit{p}}_{\Lambda |{\mathit{X}}_{1},\dots ,{\mathit{X}}_{\mathit{n}}}\left(\lambda \phantom{\rule{0.3em}{0ex}}|\phantom{\rule{0.3em}{0ex}}{\mathit{x}}_{1},\dots ,{\mathit{x}}_{\mathit{n}}\right)\text{d}\lambda \phantom{\rule{2em}{0ex}}& \hfill & \phantom{\rule{2em}{0ex}}\\ \hfill & =\frac{{\mathit{n}}^{\mathit{s}+\text{\xbd}}}{\Gamma \left(\mathit{s}+\text{\xbd}\right)}\frac{\Gamma \left(\mathit{y}+\mathit{s}+\text{\xbd}\right)}{\mathit{y}!{\left(\mathit{n}+1\right)}^{\mathit{y}+\mathit{s}+\text{\xbd}}}.\phantom{\rule{2em}{0ex}}& \hfill \text{(5)}\end{array}$$This defines a discrete probability distribution on the non-negative integers, often called a Poisson-gamma mixture distribution [Bernardo and Smith, 2000, §3.2.2]. For our data, since $\pi $ achieves a maximum at $\mathit{y}=5$ (Figure 5), this is the (a posteriori) most likely number $\mathit{Y}$ of major hurricanes that will hit the U.S. mainland in 2011–2020. The mean of the posterior distribution is $6$. Since the probability is 0.956 that $\mathit{Y}$’s value lies between 2 and 11 (inclusive), the interval whose end-points are 2 and 11 is a $95.6\phantom{\rule{0.3em}{0ex}}\%$ coverage interval for $\mathit{Y}$.

[Altman, 1991, Table 9.6] lists measurement results from Hulse et al. [1979], for the concentration of thyroxine in the serum of sixteen children diagnosed with hypothyroidism, of which nine had slight or no symptoms, and the other seven had marked symptoms. The values measured for the former, all in units of nmol/L, were 34, 45, 49, 55, 58, 59, 60, 62, and 86; and for the latter they were 5, 8, 18, 24, 60, 84, and 96. The averages are $\overline{\mathit{x}}=56.4$ and $\overline{\mathit{y}}=42.1$, and the standard deviations are $\mathit{s}=14.2$ and $\mathit{t}=37.5$.

Our goal is to produce a probability interval for the difference between the corresponding means, $\mu $ and $\nu $, say, when nothing is assumed known a priori either about these means or about the corresponding standard deviations, $\sigma $ and $\tau $, which may be different.

Given the values of these four parameters, suppose that the values measured in the $\mathit{m}=9$ children with slight or no symptoms are observed values of independent Gaussian random variables ${\mathit{U}}_{1},\dots ,{\mathit{U}}_{\mathit{m}}$ with common mean $\mu $ and standard deviation $\sigma $, and that those measured in the $\mathit{n}=7$ children with marked symptoms are observed values of independent Gaussian random variables ${\mathit{V}}_{1},\dots ,{\mathit{V}}_{\mathit{n}}$, also independent of the $\left\{{\mathit{U}}_{\mathit{i}}\right\}$, with common mean $\nu $ and standard deviation $\tau $.

The problem of constructing a probability interval for $\mu -\nu $ under these circumstances is known as the Behrens-Fisher problem [Ghosh and Kim, 2001]. For the Bayesian solution, we regard $\mu $, $\nu $, $\sigma $ and $\tau $ as realized values of non-observable random variables $\mathit{M}$, $\mathit{N}$, $\Sigma $, and $\mathit{T}$, assumed independent a priori and such that $\mathit{M}$, $\mathit{N}$, $log\Sigma $, and $log\mathit{T}$ all are uniformly distributed over the real numbers (hence have improper prior distributions). The corresponding posterior distributions all are proper provided $\mathit{m}\ge 2$ and $\mathit{n}\ge 2$. However, the density of the posterior probability distribution of $\mathit{M}-\mathit{N}$ given the data cannot be computed in closed form.

This problem in Bayesian inference, and other problems much more demanding than this, can be solved using the MCMC sampling technique mentioned in §6.4, for which there exist several generic software implementations: we obtained the results presented below using function metrop of the R package mcmc [Geyer, 2010]. Typically, all that is needed is the logarithm of the numerator of Bayes formula (3). Leaving out constants that do not involve $\mu $, $\nu $, $\sigma $ or $\tau $, this is

$$-\left(\mathit{m}+1\right)log\left(\sigma \right)-\left(\mathit{n}+1\right)log\left(\tau \right)-\frac{\mathit{m}{\left(\mu -\overline{\mathit{x}}\right)}^{2}+\left(\mathit{m}-1\right){\mathit{s}}^{2}}{2{\sigma}^{2}}-\frac{\mathit{n}{\left(\nu -\overline{\mathit{y}}\right)}^{2}+\left(\mathit{n}-1\right){\mathit{t}}^{2}}{2{\tau}^{2}}.$$ |

MCMC produces a sample of suitably large size $\mathit{K}$ from the joint posterior distribution of $\mathit{M}$, $\mathit{N}$, $\Sigma $, and $\mathit{T}$, given the data, say $\left({\mu}_{1},{\nu}_{1},{\sigma}_{1},{\tau}_{1}\right)$, …, $\left({\mu}_{\mathit{K}},{\nu}_{\mathit{K}},{\sigma}_{\mathit{K}},{\tau}_{\mathit{K}}\right)$. The 95 % probability interval for the difference in mean levels of thyroxine in the serum of the two groups, which extends from $\text{-22}\phantom{\rule{0.3em}{0ex}}\text{nmol/L}$ to $\text{51}\phantom{\rule{0.3em}{0ex}}\text{nmol/L}$, and Figure 6, are based on a sample of size $\mathit{K}=4.5\times 1{0}^{6}$. The probability density in this figure, and that probability interval, were computed as described in MC4.c and MC4.d of §5.8, only applied to the differences $\left\{{\mu}_{\mathit{k}}-{\nu}_{\mathit{k}}\right\}$.

In this particular case it is possible to ascertain the correctness of the results owing to an interesting, albeit surprising result: the posterior means are independent a posteriori, and have probability distributions that are re-scaled, shifted versions of Student’s $\mathit{t}$ distributions with $\mathit{m}-1$ and $\mathit{n}-1$ degrees of freedom [Box and Tiao, 1973, 2.5.2]. Therefore, by application of the Monte Carlo method of §5.8, one may obtain a sample from the posterior distribution of the difference $\mathit{M}-\mathit{T}$ independently of the MCMC procedure: the results are depicted in Figure 6 (where they are labeled “Jeffreys (Exact)”), and are essentially indistinguishable from the results of MCMC.

This same figure shows yet another posterior density that differs hardly at all from the posterior density corresponding to Jeffreys’s prior: this alternative result corresponds to the “matching” prior distribution (also improper) derived by Ghosh and Kim [2001], whose density is proportional to $\left({\sigma}^{2}\u2215\mathit{m}+{\tau}^{2}\u2215\mathit{n}\right)\u2215{\left(\sigma \tau \right)}^{3}$. This illustrates a generally good practice: that the sensitivity of the results of Bayesian analysis should be evaluated by comparing how they vary when different but comparably acceptable priors are used.

Antonio Possolo thanks colleagues and former colleagues from the Working Group 1 of the Joint Committee for Guides in Metrology, who offered valuable comments and suggestions in multiple discussions of an early draft of a document that he had submitted to the consideration of this Working Group and that included material now in the present document: Walter Bich, Maurice Cox, René Dybkær, Charles Ehrlich, Clemens Elster, Tyler Estler, Brynn Hibbert, Hidetaka Imai, Willem Kool, Lars Nielsen, Leslie Pendrill, Lorenzo Peretto, Steve Sidney, Adriaan van der Veen, Graham White, and Wolfgang Wöger. Antonio Possolo is particularly grateful to Tyler Estler for suggestions regarding §2.3 and §3.1, and to Graham White for suggestions and corrections that greatly improved §3.6, on the tuberculin test. Both authors thank their common NIST colleagues Andrew Rukhin and Jack Wang for suggesting many corrections and improvements. Mary Dal-Favero and Alan Heckert, both from NIST, kindly facilitated the deployment of this material on the World Wide Web.

D. G. Altman. Practical Statistics for Medical Research. Chapman & Hall/CRC, Boca Raton, FL, 1991. Reprinted 1997.

American Thoracic Society. Diagnostic standards and classification of tuberculosis in adults and children. American Journal of Respiratory and Critical Care Medicine, 161:1376–1395, 1999.

J. Berger. The case for objective Bayesian analysis. Bayesian Analysis, 1(3):385–402, 2006. URL http://ba.stat.cmu.edu/.

J. Bernardo and A. Smith. Bayesian Theory. John Wiley & Sons, New York, 2000.

J. Bernardo and A. Smith. Bayesian Theory. John Wiley & Sons, Chichester, England, 2nd edition, 2007.

J. M. Bernardo. Reference posterior distributions for Bayesian inference. Journal of the Royal Statistical Society, 41:113–128, 1979.

E. S. Blake, C. W. Landsea, and E. J. Gibney. The deadliest, costliest, and most intense United States tropical cyclones from 1851 to 2010 (and other frequently requested hurricane facts). Technical Report Technical Memorandum NWS NHC-6, NOAA, National Weather Service, National Hurricane Center, Miami, Florida, August 2011.

L. Bovens and W. Rabinowicz. Democratic answers to complex questions — an epistemic perspective. Synthese, 150:131–153, 2006.

G. E. P. Box and G. C. Tiao. Bayesian Inference in Statistical Analysis. Addison-Wesley, Reading, Massachusetts, 1973.

T. W. Cannon. Light and radiation. In C. DeCusatis, editor, Handbook of Applied Photometry, chapter 1, pages 1–32. Springer Verlag, New York, New York, 1998.

R. Carnap. Logical Foundations of Probability. University of Chicago Press, Chicago, Illinois, 2nd edition, 1962.

G. Casella and R. L. Berger. Statistical Inference. Duxbury, Pacific Grove, California, 2nd edition, 2002.

N. Clee. Who’s the daddy of them all? In Observer Sport Monthly. Guardian News and Media Limited, Manchester, UK, Sunday March 4, 2007.

R. T. Clemen and R. L. Winkler. Combining probability distributions from experts in risk analysis. Risk Analysis, 19:187–203, 1999.

R. T. Cox. Probability, frequency and reasonable expectation. American Journal of Physics, 14:1–13, 1946.

R. T. Cox. The Algebra of Probable Inference. The Johns Hopkins Press, Baltimore, Maryland, 1961.

A. C. Davison and D. Hinkley. Bootstrap Methods and their Applications. Cambridge University Press, New York, NY, 1997.

B. de Finetti. La prévision: ses lois logiques, ses sources subjectives. Anales de l’Institut Henri Poincaré, 7:1–68, 1937.

B. de Finetti. Theory of Probability: A critical introductory treatment. John Wiley & Sons, Chichester, 1990. Two volumes, translated from the Italian and with a preface by Antonio Machì and Adrian Smith, with a foreword by D. V. Lindley, Reprint of the 1975 translation.

M. H. DeGroot. A conversation with Persi Diaconis. Statistical Science, 1(3):319–334, August 1986.

M. H. DeGroot and M. J. Schervish. Probability and Statistics. Addison-Wesley, 4th edition, 2011.

P. Diaconis and D. Ylvisaker. Quantifying prior opinion. In J. M. Bernardo, M. H. DeGroot, D. V. Lindley, and A. Smith, editors, Bayesian Statistics, volume 2, pages 163–175. North-Holland, Amsterdam, 1985.

F. W. Dyson, A. S. Eddington, and C. Davidson. A determination of the deflection of light by the sun’s gravitational field, from observations made at the total eclipse of May 29, 1919. Philosophical Transactions of the Royal Society of London, Series A, Containing Papers of a Mathematical or Physical Character, 220:291–333, 1920.

B. Efron and R. J. Tibshirani. An Introduction to the Bootstrap. Chapman & Hall, London, UK, 1993.

W. Feller. An Introduction to Probability Theory and Its Applications, volume I. John Wiley & Sons, New York, 3rd edition, 1968. Revised Printing.

W. Feller. An Introduction to Probability Theory and Its Applications, volume II. John Wiley & Sons, New York, 2nd edition, 1971.

R. A. Fisher. Statistical Methods for Research Workers. Hafner Publishing Company, New York, NY, 14th edition, 1973.

B. Fitelson. Likelihoodism, bayesianism, and relational confirmation. Synthese, 156(3):473–489, 2007.

P. H. Garthwaite, J. B. Kadane, and A. O’Hagan. Statistical methods for eliciting probability distributions. Journal of the American Statistical Association, 100:680–701, June 2005.

C. F. Gauss. Theoria combinationis observationum erroribus minimis obnoxiae. In Werke, Band IV . Könighlichen Gesellschaft der Wissenschaften, Göttingen, 1823.

A. Gelman, J. B. Carlin, H. S. Stern, and D. B. Rubin. Bayesian Data Analysis. Chapman & Hall / CRC, 2nd edition, 2003.

C. J. Geyer. mcmc: Markov Chain Monte Carlo, 2010. URL http://CRAN.R-project.org/package=mcmc. R package version 0.8.

M. Ghosh and Y.-H. Kim. The Behrens-Fisher problem revisited: A Bayes-Frequentist synthesis. The Canadian Journal of Statistics, 29(1): 5–17, March 2001.

D. Gillies. Philosophical Theories of Probability. Routledge, London, UK, 2000.

C. Glymour. Theory and evidence. Princeton University Press, Princeton, New Jersey, 1980.

A. Hájek. Interpretations of probability. In E. N. Zalta, editor, The Stanford Encyclopedia of Philosophy. The Metaphysics Research Lab, Center for the Study of Language and Information, Stanford University, Stanford, California, 2007. URL http://plato.stanford.edu/archives/win2007/entries/probability-interpret/.

S. Hartmann and J. Sprenger. Judgment aggregation and the problem of tracking the truth. Synthese, pages 1–13, 2011. URL http://dx.doi.org/10.1007/s11229-011-0031-5.

P.G. Hoel, S. C. Port, and C. J. Stone. Introduction to Probability Theory. Houghton Mifflin, 1971a.

P.G. Hoel, S. C. Port, and C. J. Stone. Introduction to Statistical Theory. Houghton Mifflin, 1971b.

M. Holden, M. R. Dubin, and P. H. Diamond. Frequency of negative intermediate-strength tuberculin sensitivity in patients with active tuberculosis. New England Journal of Medicine, 285:1506–1509, 1971.

J. A. Hulse, D. Jackson, D. B. Grant, P. G. H. Byfield, and R. Hoffenberg. Different measurements of thyroid function in hypothyroid infants diagnosed by screening. Acta Pædiatrica, 68:21–25, 1979.

E. T. Jaynes. Probability Theory in Science and Engineering. Colloquium Lectures in Pure and Applied Science, No. 4. Socony Mobil Oil Company, Dallas, Texas, 1958.

E. T. Jaynes. Probability Theory: The Logic of Science. Cambridge University Press, Cambridge, UK, 2003. G. L. Bretthorst, Editor.

H. Jeffreys. Theory of Probability. Oxford University Press, London, 3rd edition, 1961. Corrected Impression, 1967.

Joint Committee for Guides in Metrology. Evaluation of measurement data — Guide to the expression of uncertainty in measurement. International Bureau of Weights and Measures (BIPM), Sèvres, France, September 2008a. URL http://www.bipm.org/en/publications/guides/gum.html. BIPM, IEC, IFCC, ILAC, ISO, IUPAC, IUPAP and OIML, JCGM 100:2008, GUM 1995 with minor corrections.

Joint Committee for Guides in Metrology. Evaluation of measurement data — Supplement 1 to the “Guide to the expression of uncertainty in measurement” — Propagation of distributions using a Monte Carlo method. International Bureau of Weights and Measures (BIPM), Sèvres, France, 2008b. URL http://www.bipm.org/en/publications/guides/gum.html. BIPM, IEC, IFCC, ILAC, ISO, IUPAC, IUPAP and OIML, JCGM 101:2008.

Joint Committee for Guides in Metrology. International vocabulary of metrology — Basic and general concepts and associated terms (VIM). International Bureau of Weights and Measures (BIPM), Sèvres, France, 2008c. URL http://www.bipm.org/en/publications/guides/vim.html. BIPM, IEC, IFCC, ILAC, ISO, IUPAC, IUPAP and OIML, JCGM 200:2008.

R. Köhler. Photometric and radiometric quantities. In C. DeCusatis, editor, Handbook of Applied Photometry, chapter 2, pages 33–54. Springer Verlag, New York, New York, 1998.

D. Lindley. Understanding Uncertainty. John Wiley & Sons, Hoboken, New Jersey, 2006.

D. V. Lindley. Introduction to Probability and Statistics from a Bayesian Viewpoint — Part 1, Probability. Cambridge University Press, Cambridge, UK, 1965a.

D. V. Lindley. Introduction to Probability and Statistics from a Bayesian Viewpoint — Part 2, Inference. Cambridge University Press, Cambridge, UK, 1965b.

D. V. Lindley. Reconciliation of probability distributions. Operations Research, 31(5):866–880, September-October 1983.

D. V. Lindley. Making Decisions. John Wiley & Sons, London, 2nd edition, 1985.

D. H. Mellor. Probability: A Philosophical Introduction. Routledge, New York, 2005.

N. Metropolis and S. Ulam. The Monte Carlo Method. Journal of the American Statistical Association, 44:335–341, September 1949.

N. Metropolis, A. Rosenbluth, M. Rosenbluth, A. Teller, and E. Teller. Equations of state calculations by fast computing machines. Journal of Chemical Physics, 21:1087–1091, 1953.

M. G. Morgan and M. Henrion. Uncertainty — A Guide to Dealing with Uncertainty in Quantitative Risk and Policy Analysis. Cambridge University Press, New York, NY, first paperback edition, 1992. 10th printing, 2007.

P. A. Morris. Combining expert judgments: A bayesian approach. Management Science, 23(7):679–693, March 1977.

J. Neyman. Frequentist probability and frequentist statistics. Synthese, 36(1):97–131, 1977.

K. R. Popper. The propensity interpretation of probability. British Journal of the Philosophy of Science, 10:25–42, 1959.

A. Possolo. Copulas for uncertainty analysis. Metrologia, 47:262–271, 2010.

R Development Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2010. URL http://www.R-project.org. ISBN 3-900051-07-0.

F. P. Ramsey. Truth and probability. In R.B. Braithwaite, editor, The Foundations of Mathematics and other Logical Essays, chapter VII, pages 156–198. Harcourt, Brace and Company, New York, 1999 electronic edition, 1926, 1931. URL http://homepage.newschool.edu/het/texts/ramsey/ramsess.pdf.

H. Reichenbach. The Theory of Probability: An Inquiry into the Logical and Mathematical Foundations of the Calculus of Probability. University of California Press, Berkeley, California, 1949. English translation of the 1935 German edition.

C. P. Robert. The Bayesian Choice: From Decision-Theoretic Foundations to Computational Implementation. Springer, New York, NY, second edition, 2007.

L. J. Savage. The Foundations of Statistics. Dover Publications, New York, New York, 1972.

M. J. Schervish. Theory of Statistics. Springer Series in Statistics. Springer Verlag, New York, NY, 1995.

B. W. Silverman. Density Estimation. Chapman and Hall, London, 1986.

R. H. Simpson. The hurricane disaster potential scale. Weatherwise, 27: 169–186, 1974.

M. Stone. The opinion pool. The Annals of Mathematical Statistics, 32: 1339–1342, December 1961.

Student. The probable error of a mean. Biometrika, 6(1):1–25, March 1908.

B. N. Taylor and C. E. Kuyatt. Guidelines for Evaluating and Expressing the Uncertainty of NIST Measurement Results. National Institute of Standards and Technology, Gaithersburg, MD, 1994. URL http://physics.nist.gov/Pubs/guidelines/TN1297/tn1297s.pdf. NIST Technical Note 1297.

R. von Mises. Probability, Statistics and Truth. Dover Publications, New York, 2nd revised edition, 1981. ISBN 0486242145. Translation of the 3rd German edition.