Chapter 2

Probability and Statistics Review

In this chapter we briefly review without proofs some definitions and concepts in probability and statistics. Many introductory and more advanced texts can be recommended for review and reference. On introductory probability see e.g. Bean [23], Ghahramani [118], or Ross [232]. Mathematical statistics and probability books at an advanced undergraduate or first year graduate level include e.g. DeGroot and Schervish [64], Freund (Miller and Miller) [201], Hogg, McKean and Craig [146] or Larsen and Marx [170]. Casella and Berger [39] or Bain and Englehart [16] are somewhat more advanced. Durrett [77] is a graduate probability text. Lehmann [172] and Lehmann and Casella [173] are graduate texts in statistical inference.

2.1 Random Variables and Probability

Distribution and Density Functions

The cumulative distribution function (cdf) of a random variable X is FX defined by

$\begin{array}{cc}{F}_X(x)=P(X\le x),& x\end{array}\in ℝ.$

In this book P(·) denotes the probability of its argument. We will omit the subscript X and write F(x) if it is clear in context. The cdf has the following properties:

1. FX is non-decreasing.
2. FX is right-continuous; that is,

   $\begin{array}{cc}\underset{\in \to 0^{+}}{\lim }{F}_x(x+\in )=F_X(x),& \text{for}\text{all}x\in ℝ.\end{array}$

3. $\begin{array}{ccc}\underset{x\to -\infty }{\lim }{F}_X(x)=0& and& \underset{x\to \infty }{\lim }{F}_X(x)=1.\end{array}$

A random variable X is continuous if FX is a continuous function. A random variable X is discrete if FX is a step function.

Discrete distributions can be specified by the probability mass function (pmf) pX(x) = P(X = x). The discontinuities in the cdf are at the points where the pmf is positive, and $p\left(x\right)=F\left(x\right)-F\left(x^{-}\right).$

If X is discrete, the cdf of X is

$F_X(x)=P(X\le x)={\displaystyle \sum _{\left\{k\le x:p_X(k)>0\right\}}{p}_X(k).}$

Continuous distributions do not have positive probability mass at any single point. For continuous random variables X the probability density function (pdf) or density of X is $f_X\left(x\right)=F_X^{\prime }\left(x\right)$, provided that FX is differentiable, and by the fundamental theorem of calculus

$F_X\text{(}x\text{) = }P\text{(}X\le x\text{) =}{\displaystyle {\int }_{-\infty }^x{f}_x(t)dt.}$

The joint density of continuous random variables X and Y is fX,Y (x, y) and the cdf of (X, Y) is

${\text{F}}_{X,Y}\text{(}x\text{, }y\text{) = }P\text{(}X\le x\text{; }Y\le y\text{) =}{{\displaystyle \int }}_{-\infty }^y{{\displaystyle \int }}_{-\infty }^x{\text{f}}_{X,Y}\text{(}s\text{, }t\text{)}dsdt.$

The marginal probability densities of X and Y are given by

$\begin{array}{cc}{f}_X(x)={{\displaystyle \int }}_{-\infty }^{\infty }{f}_{X,Y}(x,y)dy;& {\text{f}}_Y\text{ (}y\text{)}={{\displaystyle \int }}_{-\infty }^{\infty }{\text{f}}_{X,Y}\text{(}x\text{, }y\text{)}dx.\end{array}$

The corresponding formulas for discrete random variables are similar, with sums replacing the integrals. In the remainder of this chapter, for simplicity fX(x) denotes either the pdf (if X is continuous) or the pmf (if X is discrete) of X.

The set of points {x : fX(x) > 0} is the support set of the random variable X. Similarly, the bivariate distribution of (X, Y) is supported on the set {(x, y) : fX,Y (x, y) > 0}.

Expectation, Variance, and Moments

The mean of a random variable X is the expected value or mathematical expectation of the variable, denoted E[X]. If X is continuous with density f, then the expected value of X is

$E[X]={\displaystyle {\int }_{-\infty }^{\infty }xf(x)dx.}$

If X is discrete with pmf f(x), then

$E[X]={\displaystyle \sum }_{\text{{}x\text{: }{f}_X\text{(}x\text{)>0}}}xf\text{(}x\text{)}\text{.}$

(The integrals and sums above are not necessarily finite. We implicitly assume that E|X|< ∞ whenever E[X] appears in formulas below.) The expected value of a function g(X) of a continuous random variable X with pdf f is defined by

$E\text{[}g\text{(}X\text{)]}={\displaystyle {\int }_{-\infty }^{\infty }g\text{(}x\text{)}f\text{(}x\text{)}dx.}$

Let μX = E[X].Then μX is also called the first moment of X. The rth moment of X is $E[X^r]$. Hence ifX is continuous,

$E[X^r]={\displaystyle {\int }_{-\infty }^{\infty }{x}^r{f}_X(x)dx}.$

The variance of X is the second central moment,

Var(X) = E[(X − E[X])2].

The identity E[(X − E[X])2] = E[X2] − (E[X])2 provides an equivalent formula for variance,

$\begin{array}{ccc}\mathrm{var}\left(X\right)=& E\left[X^2\right]-{\left(E\left[X\right]\right)}^2=& E\left[X^2\right]-{\mu }_X^2\end{array}.$

The variance of X is also denoted by ${\sigma }_X^2$. The square root of the variance is the standard deviation. The reciprocal of the variance is the precision.

The expected value of the product of continuous random variables X and Y with joint pdf fX,Y is

$E[XY]={\displaystyle {\int }_{-\infty }^{\infty }{\displaystyle {\int }_{-\infty }^{\infty }xy{f}_{X,Y}(x,y)}}dxdy.$

The covariance of X and Y is defined by

$\begin{array}{lll}Cov\left(X,Y\right)\hfill & =\hfill & E\left[\left(X-{\mu }_X\right)\left(Y-{\mu }_Y\right)\right]\hfill \\ \hfill & =\hfill & E\left[XY\right]-E\left[X\right]E\left[Y\right]=E\left[XY\right]-{\mu }_X{\mu }_Y.\hfill \end{array}$

The covariance of X and Y is also denoted by σXY . Note that Cov(X, X) = Var(X). The product-moment correlation is

$\rho (X,Y)=\frac{Cov\text{(}X\text{, }Y\text{ )}}{\sqrt{Var(X)Var(Y)}}=\frac{{\sigma }_{XY}}{{\sigma }_X{\sigma }_Y}.$

Correlation can also be written as

$\rho \left(X,Y\right)=E\left[\left(\frac{X-{\mu }_X}{{\sigma }_X}\right)\left(\frac{Y-{\mu }_Y}{{\sigma }_Y}\right)\right].$

Two variables X and Y are uncorrelated if ρ(X,Y) = 0.

Conditional Probability and Independence

In classical probability, the conditional probability of an event A given that event B has occurred is

$P\left(A|B\right)=\frac{P\left(AB\right)}{P\left(B\right)},$

where AB = A ∩ B is the intersection of events A and B. Events A and B are independent if P(AB) = P(A)P(B); otherwise they are dependent. The joint probability that both A and B occur can be written

$P\left(AB\right)=P\left(A|B\right)P\left(B\right)=P\left(B|A\right)P\left(A\right).$

If random variables X and Y have joint density fX,Y (x, y), then the conditional density of X given Y = y is

$f_{X|Y=y^{(x)}}=\frac{f_{X,Y}(x,y)}{f_Y(y)}.$

Similarly the conditional density of Y given X = x is

$f_{Y|X=x^{\left(y\right)}}=\frac{f_{X,Y}(x,y)}{f_X(x)}.$

Thus, the joint density of (X, Y) can be written

$f_{\text{x,y}}\left(x,y\right)=f_{X\text{|}Y=y}\left(x\right)f_{\text{y}}\left(y\right)=f_{Y\text{|}X=x}\left(y\right)f_{\text{x}}\left(x\right).$

Independence

The random variables X and Y are independent if and only if

$f_{\text{X,Y}}\left(x,y\right)=f_{\text{X}}\left(x\right)f_{\text{Y}}\left(y\right)$

for all x and y; or equivalently, if and only if $F_{\text{X,Y}}\left(x,y\right)=F_{\text{X}}\left(x\right)F_{\text{Y}}\left(y\right)$, for all x and y.

The random variables X1, ... , Xd are independent if and only if the joint pdf f of X1, ... , Xd is equal to the product of the marginal density functions. That is, X1, ... , Xd are independent if and only if

$f\left(x_{1,\mathrm{...},}{x}_d\right)={\displaystyle \prod _{j=1}^d{f}_j\left(x_j\right)}$

for all x = (x1, ... , xd)T in ℝd, where fj(xj) is the marginal density (or marginal pmf) of Xj.

The variables {X1, ... ,Xn}are a random sample from a distribution FX if X1, ... , Xn are independently and identically distributed with distribution FX. In this case the joint density of {X1, ... , Xn} is

$f\left(x_{1,\mathrm{...},}{x}_n\right)={\displaystyle \prod _{i=1}^n{f}_X\left(x_i\right).}$

If X and Y are independent, then Cov(X,Y) = 0 and p(X,Y) = 0. However, the converse is not true; uncorrelated variables are not necessarily independent. The converse is true in an importantspecial case: if X and Y are normally distributed then Cov(X,Y) = 0 implies independence.

Properties of Expected Value and Variance

Suppose that X and Y are random variables, and a and b are constants. Then the following properties hold (provided the moments exist).

1. E[aX + b] = aE[X] + b.
2. E[X + Y] = E[X] + E[Y].
3. If X and Y are independent, E[XY] = E[X]E[Y].
4. Var(b) = 0.
5. Var[aX + b] = a2Var(X).
6. Var(X + Y) = Var(X) +Var(Y) + 2Cov(X,Y).
7. If X and Y are independent, Var(X + Y) = Var(X) +Var(Y).

If {X1, ... , Xn} are independent and identically distributed (iid) we have

E[X1 + · · · + Xn] = nμX, Var(X1 + · · · + Xn) = nσ2X,

so the sample mean $\overline{X}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{X}_i}$ has expected value μX and variance $\sigma {}_X^2/n$. (Apply properties 2, 7, and 5 above.)

The conditional expected value of X given Y = y is

$E\left[X|Y=y\right]={\displaystyle {\int }_{-\infty }^{\infty }x{f}_{X|Y=y}\left(x\right)dx,}$

if $f_{X|Y=y}\left(x\right)$ is continuous.

Two important results are the conditional expectation rule and the conditional variance formula:

$\begin{array}{cc}E\left[X\right]=E\left[E\left[X|Y\right]\right]& \left(2.1\right)\end{array}$

$\begin{array}{cc}\mathrm{var}\left(X\right)=E\left[\mathrm{var}\left(X|Y\right)\right]+\mathrm{var}\left(E\left[X|Y\right]\right)& \left(2.2\right)\end{array}$

See e.g. Ross [233, Ch. 3] for a proof of (2.1, 2.2) and many applications.

2.2 Some Discrete Distributions

Some important discrete distributions are the “counting distributions.” The counting distributions are used to model the frequency of events and waiting time for events in discrete time, for example. Three important counting distributions are the binomial (and Bernoulli), negative binomial (and geometric), and Poisson.

Several discrete distributions including the binomial, geometric, and negative binomial distributions can be formulated in terms of the outcomes of Bernoulli trials. A Bernoulli experiment has exactly two possible outcomes, “success” or “failure.” A Bernoulli random variable X has the probability mass function

P(X = 1) = p, P(X = 0) = 1 − p,

where p is the probability of success. It is easy to check that E[X] = p and Var(X) = p(1 − p). A sequence of Bernoulli trials is a sequence of outcomes X1, X2,... of iid Bernoulli experiments.

Binomial and Multinomial Distribution

Suppose that X records the number of successes in n iid Bernoulli trials with success probability p. Then X has the Binomial(n, p) distribution [abbreviated X ~ Bin(n, p)] with

$P\left(X=x\right)=\left(\begin{array}{c}n\\ x\end{array}\right)p^x{\left(1-p\right)}^{n-x}=\frac{n!}{x!\left(n-x\right)!}p{\left(1-p\right)}^{n-x},x=0,1,\mathrm{...},n.$

The mean and variance formulas are easily derived by observing that that the binomial variable is an iid sum of n Bernoulli(p) variables. Therefore

E[X] = np,       Var(X) = np(1 − p).

A binomial distribution is a special case of a multinomial distribution. Suppose that there are k+1 mutually exclusive and exhaustive events A1,... , Ak+1 that can occur on any trial of an experiment, and each event occurs with probability P(Aj) = pj, j = 1, ... , k + 1. Let Xj record the number of times that event Aj occurs in n independent and identical trials of the experiment. Then X = (X1,..., Xk) has the multinomial distribution with joint pdf

$\begin{array}{l}f\left(x_{1,\mathrm{...},}{x}_k\right)=\frac{n!}{x_1!x_2!\mathrm{...}{x}_{k+1}!}{p}_1^{x_1}{p}_2^{x_2}\mathrm{...}{p}_{k+1}^{x_{k+1}}0\le x_j\le n,\left(2.3\right)\\ \text{where}{x}_{k+1}=n-{\displaystyle {\sum }_{j=1}^k{x}_j}.\end{array}$

Geometric Distribution

Consider a sequence of Bernoulli trials, with success probability p. Let the random variable X record the number of failures until the first success is observed. Then

$\begin{array}{cc}P\left(X=x\right)=p{\left(1-p\right)}^x,x=0,1,2....& \left(2.4\right)\end{array}$

A random variable X with pmf (2.4) has the Geometric(p) distribution [abbreviated X ~ Geom(p)]. If X ~ Geom(p), then the cdf of X is

$F_X\left(x\right)=P\left(X\le x\right)=1-{\left(1-p\right)}^{\left[x\right]+1},x\ge 0,$

and otherwise FX(x) = 0. The mean and variance of X are given by

$E\left[X\right]=\frac{1-p}{p};\text{V ar[}X\text{]=}\frac{1-p}{p^2}.$

Alternative formulation of Geometric distribution

The geometric distribution is sometimes formulated by letting Y be defined as the number of trials until the first success. Then Y = X + 1, where X is the random variable defined above with pmf (2.4). Under this model, we have P(Y = y) = p(1 − p)y−1, y = 1, 2, ... , and

$E\text{[}Y\text{ ] = }E\text{[}X\text{ + 1] =}\frac{1-p}{p}+1=\frac{1}{p};$

$\text{V ar[}Y\text{ ] = V ar[}X\text{ + 1] = V ar[}X\text{] =}\frac{1-p}{p^2}.$

However, as a counting distribution, or frequency model, the first formulation (2.4) given above is usually applied, because frequency models typically must include the possibility of a zero count.

Negative Binomial Distribution

The negative binomial frequency model applies in the same setting as a geometric model, except that the variable of interest is the number of failures until the r success. Suppose that exactly X failures occur before the rth success. If X = x, then the r success occurs on the (x + r)th trial. In the first x + r − 1 trials, there are r − 1 successes and x failures. This can happen $\left(\begin{array}{c}x+r-1\\ r-1\end{array}\right)=\left(\begin{array}{c}x+r-1\\ x\end{array}\right)$ ways, and each way has probability prqx. The probability mass function of the random variable X is given by

$P\left(X=x\right)=\left(\begin{array}{c}x+r-1\\ r-1\end{array}\right)p^r{q}^x,x=0,1,2,\mathrm{....}\left(2.5\right)$

The negative binomial distribution is defined for r > 0 and 0 < p < 1 as follows. The random variable X has a negative binomial distribution with parameters (r, p) if

$P\left(X=x\right)=\frac{\Gamma \left(x+r\right)}{\Gamma \left(r\right)\Gamma \left(x+1\right)}{p}^r{q}^x,x=0,1,2,\mathrm{...},\left(2.6\right)$

where Γ(⋅) is the complete gamma function defined in (2.8). Note that (2.5)and (2.6) are equivalent when r is a positive integer. If X has pmf (2.6) we will

write X ~ NegBin(r,p). The special case NegBin(r = 1,p) is the Geom(p) distribution.

Suppose that X ~ NegBin(r, p), where r is a positive integer. Then X is the iid sum of r Geom(p) variables. Therefore, the mean and variance of X given by

$E\left[X\right]=r\frac{1-p}{p}Var\text{[}X\text{] = }r\frac{1-p}{p^2},$

are simply r times the mean and variance of the Geom(p) variable in (2.4). These formulas are also valid for all r > 0.

Note that like the geometric random variable, there is an alternative formulation of the negative binomial model that counts the number of trials until the rth success.

Poisson Distribution

A random variable X has a Poisson distribution with parameter λ > 0 if the pmf of X is

$p\left(x\right)=\frac{e^{-\text{λ}}{\text{λ}}^x}{x!},x=0,1,\mathrm{2....}$

If X ~ Poisson(λ) then

$E\left[X\right]=\text{λ;}\text{var}\left(X\right)=\text{λ}\text{.}$

A useful recursive formula for the pmf is $p\left(x+1\right)=p\left(x\right)\frac{\text{λ}}{x+1},x=0,1,2,\mathrm{....}$ The Poisson distribution has many important properties and applications (see e.g. [124, 158, 233]).

Examples

Example 2.1 (Geometric cdf)

The cdf of the geometric distribution with success probability p can be derived as follows. If q = 1− p, then at the points x = 0, 1, 2,... the cdf of X is given by

$P(X\le x)={\displaystyle \sum _{k=0}^{\infty }p{q}^k=p(1+q+q^2+\cdot \cdot \cdot +q^x)}=\frac{p(1-q^{x+1})}{1-q}=1-q^{x+1.}$

Alternately,$P(X\le x)=1-P(X\ge x+1)=1-P$(first x+1 trials are failures)$\text{=}1-q^{x+1}$.

Example 2.2 (Mean of the the Poisson distribution)

If X ~ Poisson(λ), then

$E\left[X\right]={\displaystyle \sum _{x=0}^{\infty }\frac{e^{-\text{λ}}{\text{λ}}^x}{(x)!}}=\text{λ}{\displaystyle \sum _{x=1}^{\infty }\frac{e^{-\text{λ}}{\text{λ}}^{x-1}}{(x-1)!}}=\text{λ}{\displaystyle \sum _{x=0}^{\infty }\frac{e^{-\text{λ}}{\text{λ}}^x}{(x)!}}=\text{λ}\text{.}$

The last equality follows because the summand is the Poisson pmf and the total probability must sum to 1.

2.3 Some Continuous Distributions

Normal Distribution

The normal distribution with mean μ and variance σ2 [abbreviated N(μ, σ2)] is the continuous distribution with pdf

$f(x)=\frac{1}{\sqrt{2\pi \sigma }}\exp \left\{-\frac{1}{2}{\left(\frac{x-\mu }{\sigma }\right)}^2\right\},-\infty <x<\infty .$

The standard normal distribution N(0,1) has zero mean and unit variance, and the standard normal cdf is

$\Phi (z)={\displaystyle {\int }_{-\infty }^z\frac{1}{\sqrt{2\pi }}}{e}^{-t^2/2dt,}-\infty <z<\infty .$

The normal distribution has several important properties. We summarize some of these properties, without proof. For more properties and characterizations see [156, Ch. 13], [210], or [270].

A linear transformation of a normal variable is also normally distributed. If X ~ N(μ, σ) then the distribution of Y = aX + b is N(aμ + b, a2σ2). It follows that if X ~ N(μ, σ), then

$Z=\frac{X-\mu }{\sigma }~N\left(0,1\right).$

Linear combinations of normal variables are normal; if X1, ... , Xk are independent, Xi ~ N(μi, ${\sigma }_i^2$ ), and a1, ... , ak are constants, then

Y = a1X1 + ··· + akXk

is normally distributed with mean $\mu ={\displaystyle {\sum }_{i=1}^k{a}_i{\mu }_i}$ and variance ${\sigma }^2={\displaystyle {\sum }_{i=1}^k{a}_i^2}{\sigma }_i^2.$

.

Therefore, if X1, ... , Xn is a random sample (X1, ... , Xn are iid) from a N(μ, σ2) distribution, the sum Y = X1+· · ·+Xn is normally distributed with E[Y] = nμ and Var(Y) = nσ2. It follows that the sample mean $\overline{X}=Y/n$ has the N(μ, σ2/n) distribution if the sampled distribution is normal. (In case the sampled distribution is not normal, but the sample size is large, the Central Limit Theorem implies that the distribution of Y is approximately normal. See Section 2.5)

Gamma and Exponential Distributions

A random variable X has a gamma distribution with shape parameter r > 0 and rate parameter λ > 0 if the pdf of X is

$f(x)=\frac{{\text{λ}}^r}{\Gamma (r)}{x}^{r-1}{e}^{-\text{λ}x},x\ge 0,\left(2.7\right)$

where Γ(r) is the complete gamma function, defined by

$\Gamma \left(r\right)={\displaystyle {\int }_0^{\infty }{t}^{r-1}}{e}^{-t}dt,r\ne o,-1,-2,\mathrm{....}\left(2.8\right)$

Recall thatΓ(n) = (n — 1)! for positive integers n.

The notation X ~ Gamma(r, λ) indicates that X has the density (2.7), with shape r and rate λ. If X ~ Gamma(r, λ) then

$E\left[X\right]=\frac{r}{\text{λ}};Var\text{(}X\text{) =}\frac{r}{{\text{λ}}^2}.$

Gamma distributions can also be parameterized by the scale parameter θ = 1/λ instead of the rate parameter λ. In terms of (r, θ) the mean is rθ and the variance is rθ2. An important special case of the gamma distribution is r = 1, which is the exponential distribution with rate parameter λ. The Exponential(λ) pdf is

$f\left(x\right)=\lambda e^{-\lambda x},x\ge 0$

If X is exponentially distributed with rate λ [abbreviated X ~ Exp(λ)], then

$E\left[X\right]=\frac{1}{\text{λ}};Var\text{(}X\text{) =}\frac{1}{{\text{λ}}^2}.$

It can be shown that the sum of iid exponentials has a gamma distribution.

If X1, ... , Xr are iid with the Exp(λ) distribution, then Y = X1 +... + Xr has the Gamma(r, λ) distribution.

Chisquare and t

The Chisquare distribution with ν degrees of freedom is denoted by X2(ν). The pdf of a X2(ν) random variable X is

$\text{f(x)}=\frac{1}{\Gamma \left(v/2\right)2^{(v/2)}}{x}^{(v/2)-1}{e}^{-x/2},x\in ℝ,v=1,2,\mathrm{...},.$

Note that X2(ν) is a special case of the gamma distribution, with shape parameter ν/2 and rate parameter 1/2. The square of a standard normal variable has the X2(1) distribution. If Z1,... , Zν are iid standard normal then Z21 +⋅⋅⋅+ Z2νX2(ν). If X ~X2(ν1) and Y~ X2(ν2) are independent, then X + Y ~ X2(ν1 + ν2). If X~ X2(ν), then

$E\left[X\right]=v,\mathrm{var}\left(X\right)=2v.$

The Student’s t distribution [256] is defined as follows. Let Z ~ N(0,1) and V ~ X2(ν). If Z and V are independent, then the distribution of

$T=\frac{Z}{\sqrt{V/v}}$

has the Student’s t distribution with ν degrees of freedom, denoted t(ν). The density of a t(ν) random variable X is given by

$\begin{array}{cc}f(x)=\frac{\Gamma (\frac{v+1}{2})}{\Gamma \left(\frac{v}{2}\right)}\frac{1}{\sqrt{v\pi }}\frac{1}{{\left(1+\frac{x^2}{v}\right)}^{(v+1)/2}},& x\in ℝ,v=1,2,\mathrm{...}\end{array}$

The mean and variance of X ~ t(ν) are given by

$\begin{array}{cccc}E\left[X\right]=0,& v>1;& Var(X)=\frac{v}{v-2},& v>2.\end{array}$

In the special case ν = 1 the t(1) distribution is the standard Cauchy distribution. For small ν the t distribution has “heavy tails” compared to the normal distribution. For large ν, the t(ν) distribution is approximately normal, and t(ν) converges in distribution to standard normal as ν → ∞.

Beta and Uniform Distributions

A random variable X with density function

$\begin{array}{ccc}f(x)=\frac{\Gamma \left(\alpha +\beta \right)}{\Gamma \left(\alpha \right)\Gamma \left(\beta \right)}{x}^{\alpha -1}{\left(1-x\right)}^{\beta -1},& 0\le x\le 1,\alpha >0,\beta >0.& \left(2.9\right)\end{array}$

has the Beta(α, β) distribution. The constant in the beta density is the reciprocal of the beta function, defined by

$B(\alpha ,\beta )={\displaystyle {\int }_0^1{t}^{\alpha -1}{(1-t)}^{\beta -1}}dt=\frac{\Gamma (\alpha )\Gamma \left(\beta \right)}{\Gamma \left(\alpha +\beta \right)}.$

The continuous uniform distribution on (0,1) or Uniform(0,1) is the special case Beta(1,1).

The parameters α and β are shape parameters. When α = β the distribution is symmetric about 1/2. When α ≠ β the distribution is skewed, with the direction and amount of skewness depending on the shape parameters. The mean and variance are

$\begin{array}{cc}E\left[X\right]=\frac{\alpha }{\alpha +\beta };& Var\text{(}X\text{) =}\frac{\alpha \beta }{{\left(\alpha +\beta \right)}^2\left(\alpha +\beta +1\right)}.\end{array}$

If X ~ Uniform(0, 1) = Beta(1, 1), then $E\left[X\right]=\frac{1}{2}$ and $Var\left(X\right)=\frac{1}{12}.$

In Bayesian analysis, a beta distribution is often chosen to model the distribution of a probability parameter, such as the probability of success in Bernoulli trials or a binomial experiment.

Lognormal Distribution

A random variable X has the Lognormal(μ,σ2) distribution [abbreviated X~ LogN(μ, σ2)] if X = eY , where Y ~ N(μ, σ2). That is, log X ~N(μ, σ2). The lognormal density function is

$\begin{array}{cc}{f}_X(x)=\frac{1}{x\sqrt{2\pi \sigma }}{e}^{-{(\log x-\mu )}^2/(2{\sigma }^2)},& x>0.\end{array}$

The cdf can be evaluated by the normal cdf of log X ~ N(μ,σ2), so the cdf of X ~ LogN(μ, σ2) is given by

$F_X\left(x\right)=\text{Φ}\left(\frac{\log x-\mu }{\sigma }\right),x>0$

The moments are

$\begin{array}{ccc}E\left[X^r\right]=E\left[e^{rY}\right]=\exp \left\{r\mu +\frac{1}{2}{r}^2{\sigma }^2\right\},& r>0.& (2.10)\end{array}$

The mean and variance are

$\begin{array}{cc}E\left[X\right]=e^{\mu +{\sigma }^2/2},& Var(X)=e^{2\mu +{\sigma }^2}(e^{{\sigma }^2}-1).\end{array}$

Examples

Example 2.3(Two-parameter exponential cdf)

The two-parameter exponential density is

$\begin{array}{cc}f\left(x\right)=\lambda e^{-\lambda \left(x-n\right)},x\ge n& \left(2.11\right)\end{array}$

where λ and η are positive constants. Denote the distribution with density function (2.11) by Exp(λ, η). When η = 0 the density (2.11) is exponential with rate λ.

The cdf of the two-parameter exponential distribution is given by

$\begin{array}{cc}F(x)=\underset{\eta }{\overset{x}{{\displaystyle \int }}}\text{λ}{e}^{-\text{λ(}t-\eta \text{)}}dt=\underset{0}{\overset{x-\eta }{{\displaystyle \int }}}\text{λ}{e}^{-\text{λ}u}du=1-e^{-\text{λ(}x-\eta \text{)}},& x\ge \eta .\end{array}$

In the special case η = 0 we have the cdf of the Exp(λ) distribution,

$F\left(x\right)=1-e^{-\lambda x},x\ge 0.$

Example 2.4(Memoryless property of the exponential distribution)

The exponential distribution with rate parameter λ has the memoryless property. That is, if X ~ Exp(λ), then

$P\left(X>s+t|X>s\right)=P\left(X>t\right),\text{for all }s,t\ge 0.$

The cdf of X is F(x) = 1 − exp(−λx), x ≥ 0 (see Example 2.3). Therefore, for all s, t ≥ 0 we have

$\begin{array}{l}P\left(X>s+t|X>s\right)=\frac{P\left(X>s+t\right)}{P\left(X>s\right)}=\frac{1-F\left(s+t\right)}{1-F\left(s\right)}\\ =\frac{e^{-\lambda \left(s+t\right)}}{e^{-\lambda s}}=e^{-\lambda t}=1-F\left(t\right)\\ =P\left(X>t\right).\end{array}$

The first equality is simply the definition of conditional probability, P(A|B) =P(AB)/P(B)

2.4 Multivariate Normal Distribution

The bivariate normal distribution

Two continuous random variables X and Y have a bivariate normal distribution if the joint density of (X, Y) is the bivariate normal density function, which is given by

$\begin{array}{cc}f(x,y)=\frac{1}{2\pi {\sigma }_1{\sigma }_2\sqrt{1-{\rho }^2}}\exp \left\{-\frac{1}{2\left(1-{\rho }^2\right)}\left[{\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)}^2-2p\left(\frac{x-{\mu }_1}{{\sigma }_1}\right)\left(\frac{y-{\mu }_2}{{\sigma }_2}\right)+{\left(\frac{y-{\mu }_2}{{\sigma }_2}\right)}^2\right]\right\},& (2.12)\end{array}$

(x, y) ∈ ℝ2. The parameters are μ1 = E[X], μ2 = E[Y] ${\sigma }_{{}^1}^2=\mathrm{var}\left(X\right)$,${\sigma }_{{}^2}^2=\mathrm{var}\left(Y\right)$, and ρ = Cor(X, Y). The notation (X, Y) ~ BVN(μ1, μ2, ${\sigma }_{{}^1}^2$, ${\sigma }_{{}^2}^2$,ρ) indicates that (X, Y) have the joint pdf (2.12). Some properties of the bivariate normal distribution (2.12) are:

1. The marginal distributions of X and Y are normal; that is X ~ N(μ1, ${\sigma }_{{}^1}^2$ and Y ~ N(μ2, ${\sigma }_{{}^2}^2$.
2. The conditional distribution of Y given X = x is normal with mean μ2 + ρσ2/σ1(x − μ1) and variance ${\sigma }_{{}^2}^2$ (1 − ρ2).
3. The conditional distribution of X given Y = y is normal with mean μ1 + ρσ1/σ2(y − μ2) and variance ${\sigma }_{{}^1}^2$ (1 − ρ2).
4. X and Y are independent if and only if ρ = 0.

Suppose (X1, X2) ~ BVN(μ1, μ2, ${\sigma }_{{}^1}^2$,${\sigma }_{{}^2}^2$ , ρ). Let μ = (μ1, μ2)T and

${\displaystyle \sum =\left[\begin{array}{cc}{\sigma }_{11}& {\sigma }_{12}\\ {\sigma }_{21}& {\sigma }_{22}\end{array}\right]},$

where σij = Cov(Xi,Xj). Then the bivariate normal pdf (2.12) of (X1, X2) can be written in matrix notation as rll

$f(x_1,x_2)=\frac{1}{\left(2\pi \right){\left|{\displaystyle \sum }\right|}^{1/2}}\exp \left\{-\frac{1}{2}{(x-\mu )}^T{{\displaystyle \sum \text{ }}}^{-1}(x-\mu )\right\},$

where $x={\left(x_1,x_2\right)}^T\in {ℝ}^2.$

The multivariate normal distribution

The joint distribution of continuous random variables X1, ... , Xd is multivariate normal or d-variate normal, denoted Nd(μ,$\Sigma$ ), if the joint pdf is given by

$\begin{array}{cc}f\left(x_1,...,x_d\right)=\frac{1}{{\left(2\pi \right)}^{d/2}{\left|{{\displaystyle \sum }}^{\text{ }}\right|}^{1/2}}\exp \left\{-\frac{1}{2}{\left(x-\mu \right)}^T{{{\displaystyle \sum }}^{\text{ }}}^{-1}(x-\mu )\right\},& (2.13)\end{array}$

where ∑  is the d × d nonsingular covariance matrix of (X1, ... , Xd)T, μ = (μ1, ... , μd)T is the mean vector, and x = (x1, ... , xd)T ∈ ℝd.

The one-dimensional marginal distributions of a multivariate normal variable are normal with mean μi and variance ${\sigma }_i^2$, i = 1, ... , d. Here ${\sigma }_i^2$ is the ith entry on the diagonal of ∑. In fact, all of the marginal distributions of a multivariate normal vector are multivariate normal (see e.g. Tong [273, Sec. 3.3]).

The normal random variables X1, ... , Xd are independent if and only if the covariance matrix ∑ is diagonal.

Linear transformations of multivariate normal random vectors are multivariate normal. That is, if C is an m × d matrix and b = (b1, ... , bm)T∈ ℝm, then Y = CX + b has the m-dimensional multivariate normal distribution with mean vector Cμ + b and covariance matrix C∑CT.

Applications and properties of the multivariate normal distribution are covered by Anderson [8] and Mardia et al. [188]. Refer to Tong [273] for properties and characterizations of the bivariate normal and multivariate normal distribution.

2.5 Limit Theorems

Laws of Large Numbers

The Weak Law of Large Numbers (WLLN) or (LLN) states that the sample mean converges in probability to the population mean. Suppose that X1, X2... are independent and identically distributed (iid), E|X1| < ∞ and μ = E[X1]. For each n let ${\overline{X}}_n=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{X}_i}.$ Then ${\overline{X}}_n\to \mu$ almost surely as n →∞ That is, for every ∈>0

$\underset{n\to 0}{\lim }P\left(\left|\overline{X_n}-\mu \right|<\in \right)=1.$

For a proof, see e.g. Durrett [77].

The Strong Law of Large Numbers (SLLN) states that the sample mean converges almost surely to the population mean μ. Suppose that X1, X2,... are pairwise independent and identically distributed, E|X1| < ∞ and μ = E [X1]For each n For each n let ${\overline{X}}_n=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{X}_i}.$ Then ${\overline{X}}_n\to \mu$ almost surely as n →∞ That is, for every ∈>0

$P\left(\underset{n\to 0}{\lim }\left|\overline{X_n}-\mu \right|<\in \right)=1.$

For Etemadi’s proof see Durrett [77].

Central Limit Theorem

The first version of the Central Limit Theorem was proved by de Moivre in the early 18th century for random samples of Bernoulli variables. The general proof was given independently by Lindeberg and Lévy in the early 1920’s.

THEOREM 2.1 (The Central Limit Theorem)

If X1, ... , Xn is a random sample from a distribution with mean μ and finite variance σ2> 0, then the limiting distribution of

$Z_n=\frac{\overline{X}-\mu }{\sigma /\sqrt{n}}$

is the standard normal distribution.

See Durrett [77] for the proofs.

2.6 Statistics

Unless otherwise stated, X1, ... , Xn is a random sample from a distribution with cdf FX(x) = P(X < x), pdf or pmf fX(x), mean E[X] = μX and variance ${\sigma }_X^2$. The subscript X on F, f, μ, and σ is omitted when it is clear in context. Lowercase letters x1,... , xn denote an observed random sample.

A statistic is a function Tn = T(X1, ... , Xn) of a sample. Some examples of statistics are the sample mean, sample variance, etc. The sample mean is $\overline{X}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{X}_{i,}}$ and sample variance is

$S^2=\frac{1}{n-1}{{\displaystyle \sum _{i=1}^n\left(X_i-\overline{X}\right)}}^2=\frac{{\displaystyle {\sum }_{i=1}^n{X}_i^2-n{\overline{X}}^2}}{n-1}.$

The sample standard deviation is $S=\sqrt{S^2}.$

The empirical distribution function

An estimate of F(x) = P(X ≤ x) is the proportion of sample points that fall in the interval (−∞, x]. This estimate is called the empirical cumulative distribution function (ecdf) or empirical distribution function (edf). The ecdf of an observed sample x1, ... , xn is defined by

$F_n(x)=\{\begin{array}{l}0,x<x_{(1),}\\ \frac{i}{n},x_{(i)}\le x_{(i+1),}i=1,\mathrm{...},n-1,\\ 1,x_{(n)}\le x,\end{array}$

where x(1)≤ x(2) ≤ · · · ≤ x(n) is the ordered sample.

A quantile of a distribution is found by inverting the cdf. The cdf may not be strictly increasing, however, so the definition is as follows. The q quantile of a random variable X with cdf F(x) is

$\begin{array}{cc}{X}_q=\underset{x}{\inf }\left\{x:F(x)\ge q\right\},& 0<q<1.\end{array}$

Quantiles can be estimated by the inverse ecdf of a random sample or other function of the order statistics. Methods for computing sample quantiles differ among statistical packages R, SAS, Minitab, SPSS, etc. (see Hyndman and Fan [148] and the quantile help topic in R).

R note 2.1The default method of estimation used in the R quantile function assigns cumulative probability (k − 1)/(n − 1) to the kth order statistic. Thus, the empirical cumulative probabilities are defined

$0,\frac{1}{n-1},\frac{2}{n-1},\mathrm{...},\frac{n-2}{n-1},1.$

Note that this set of probabilities differs from the usual assignment${\left\{k/n\right\}}_{k=1}^n$ of the ecdf.

Bias and Mean Squared Error

A statistic ${\widehat{\theta }}_n$ is an unbiased estimator of a parameter θ if $E[{\widehat{\theta }}_n]=\theta$ An estimatoris ${\widehat{\theta }}_n$asymptotically unbiased for θ if

$\underset{n\to \infty }{\lim }E[{\widehat{\theta }}_n]=\theta$

The bias of an estimator $\widehat{\theta }$ for a parameter θ is defined bias$(\widehat{\theta })=E[\widehat{\theta }]-\theta .$ Clearly $\overline{X}$ is an unbiased estimator of the mean $\mu =E[\widehat{\theta }]-\theta .$. It can be shown that E[S2] = σ2 = Var(X), so the sample variance S2 is an unbiased estimator of σ2. The maximum likelihood estimator of variance is

${\widehat{\sigma }}^2=\frac{1}{n}{\displaystyle \sum }_{i=1}^n{\left(X_i-\overline{X}\right)}^2$

which is a biased estimator of σ2. However, the bias −σ2/n tends to zero as n → ∞, so ${\widehat{\sigma }}^2$ is asymptotically unbiased for σ2.

The mean squared error (MSE) of an estimator $\widehat{\theta }$ for parameter θ is

$MSE\left(\widehat{\theta }\right)=E\left[{\left(\widehat{\theta }-\theta \right)}^2\right].$

Notice that for an unbiased estimator the MSE is the equal to the variance of the estimator. If $\widehat{\theta }$ is biased for θ, however, the MSE is larger than the variance. In fact, the MSE can be split into two parts,

$\begin{array}{c}MSE\left(\widehat{\theta }\right)=E\left[{\widehat{\theta }}^2-2\theta \widehat{\theta }+{\theta }^2\right]=c-2\theta E\left[\widehat{\theta }\right]+{\theta }^2\\ =E\left[{\widehat{\theta }}^2\right]-{\left(E\left[\widehat{\theta }\right]\right)}^2+{\left(E\left[\widehat{\theta }\right]\right)}^2-2\theta E\left[\widehat{\theta }\right]+{\theta }^2\\ =\mathrm{var}\left(\widehat{\theta }\right)+{\left(E\left[\widehat{\theta }\right]-\theta \right)}^2,\end{array}$

so the MSE is the sum of variance and squared bias:

$MSE(\widehat{\theta })=Var(\widehat{\theta })+{\left[bias(\widehat{\theta })\right]}^2$

The standard error of an estimator $\widehat{\theta }$ is the square root of the variance:$se(\widehat{\theta })=\sqrt{Var(\widehat{\theta })}$. An important example is the standard error of the mean

$se(\overline{X})=\sqrt{Var(\overline{X})}=\sqrt{\frac{Var(X)}{n}}=\frac{{\sigma }_x}{\sqrt{n}}$

A sample proportion $\widehat{p}$. The standard error of a sample proportion is $\sqrt{p(1-p)/n}$. Note that $se(\widehat{p})\le 0.5/\sqrt{n}$.

For each fixed x ∈ℝ, the ecdf Fn(x) is an unbiased estimator of the cdf F (x).The standard error of Fn(x) is $\sqrt{F(x)(1-F(x))/n}\le 0.5/\sqrt{n}$

The variance of the q sample quantile [63, 2.7] is

$\begin{array}{cc}Var({\widehat{x}}_q)=\frac{q(1-q)}{nf{(x_q)}^2},& (2.14)\end{array}$

where f is the density of the sampled distribution. When quantiles are estimated, the density f is usually unknown, but (2.14) shows that larger samples are needed for estimates of quantiles in the part of the support set where the density is close to zero.

Method of Moments

The rth sample moment $m_r^{’}=\frac{1}{n}{\displaystyle {\sum }_{i=1}^n{X}_i^r},$, r = 1, 2,... is an unbiased estimator of the r population moment E[X], provided that the r moment exists. If X has density f(x; θ1, ... , θk), then the method of moments estimator of θ = (θ1, ... , θk) is given by the simultaneous solution $\widehat{\theta }=\left({\widehat{\theta }}_1,...,{\widehat{\theta }}_k\right)$ of the equations

$\begin{array}{cc}E[X^r]=m_r^{’}(x_1,\mathrm{...}{x}_n)=\frac{1}{n}{\displaystyle \sum _{i=1}^n{x}_i^r},& r=1,\mathrm{...},k.\end{array}$

The Likelihood Function

Suppose that the sample observations are iid from a distribution with density function f(X|θ), where θ is a parameter. The likelihood function is the conditional probability of observing the sample, given θ, which is given by

$\begin{array}{cc}L(\theta )={\displaystyle \prod _{i=1}^{n}f(x_i|\theta ).}& (2.15)\end{array}$

The parameter θ could be a vector of parameters, θ = (θ1, ... , θp). The likelihood function regards the data as a function of the parameter(s) θ. As L(θ) is a product, it is usually easier to work with the logarithm of L(θ), called the log likelihood function,

$\begin{array}{cc}l(\theta )=\text{log}(L(\theta ))={\displaystyle \sum }_{i=1}^n\log f(x_i|\theta ).& (2.16)\end{array}$

Maximum Likelihood Estimation

The method of maximum likelihood was introduced by R. A. Fisher. By maximizing the likelihood function L(θ) with respect to θ, we are looking for the most likely value of θ given the information available, namely the sample data. Suppose that Θ is the parameter space of possible values of θ. If the maximum of L(θ) exists and it occurs at a unique point $\widehat{\theta }$∈ Θ, then $\widehat{\theta }$ is called the maximum likelihood estimator of L(θ). If the maximum exists but is not

unique, then any of the points where the maximum is attained is an MLE of θ. For many problems, the MLE can be determined analytically. However, it is often the case that the optimization cannot be solved analytically, and in that case numerical optimization or other computational approaches can be applied.

Maximum likelihood estimators have an invariance property. This property states that if $\widehat{\theta }$ is an MLE of θ and τ is a function of θ, then $\tau \left(\widehat{\theta }\right)$ is an MLE of τ(θ).

Note that the maximum likelihood principle can also be applied in problems where the observed variables are not independent or identically distributed (the likelihood function (2.15) given above is for the iid case).

Example 2.5 (Maximum likelihood estimation of two parameters)

Find the maximum likelihood estimator of θ = (λ, η) for the two-parameter exponential distribution (see Example 2.3). Suppose that x1, ... , xn is a random sample from the Exp(λ, η) distribution. The likelihood function is

$L(\theta )=L(\text{λ,}\eta )={\displaystyle \prod _{i=1}^n\text{λ}}{e}^{-\text{λ(}{x}_i\text{-}\eta \text{)}}I(x_i\ge \eta )$

where I(·) is the indicator variable (I(A) = 1 on set A and I(A) = 0 on the complement of A). Then if x(1) = min{x1, ... , xn}, we have

$\begin{array}{cc}L(\theta )=L(\text{λ,}\eta )={\text{λ}}^n\exp \left\{-\text{λ}{\displaystyle \sum _{i=1}^n(x_i-\eta )}\right\},& x_{(1)}\ge \eta ,\end{array}$

and the log-likelihood is given by

$\begin{array}{cc}l(\theta )=l(\text{λ,}\eta )=\eta \log \text{λ}-\text{λ}{\displaystyle \sum _{i=1}^n(x_i-\eta ),}& x(1)\ge \eta \end{array}$

Then l(θ) is an increasing function of η for every fixed λ, and η ≤ x(1), so $\widehat{n}=x_{\left(1\right)}$. To find the maximum of l(θ) with respect to λ, solve

$\frac{\partial l(\text{λ},\eta )}{\partial \text{λ}}=\frac{\eta }{\text{λ}}-{\displaystyle \sum _{i=1}^n(x_i-\eta )=0,}$

to find the critical point λ = 1/(x − η). The MLE of θ = (λ, η) is

$(\widehat{\lambda },\widehat{\eta })=\left(\frac{1}{\overline{x}-x_{\left(1\right)}},x(1)\right).$

Example 2.6 (Invariance property of MLE)

Find the maximum likelihood estimator of the α-quantile of the Exp(λ, η) distribution in Examples 2.3 and 2.5. From Example 2.3 we have

$F\left(x\right)=1-e^{-\lambda \left(x-n\right)},x\ge n.$

Therefore F(xα) = α implies that

$x_{\text{α}}=-\frac{1}{\text{λ}}\text{log(1-α)+}\eta ,$

and by the invariance property of maximum likelihood, the MLE of xa is

${\widehat{x}}_{\text{α}}=-(\overline{x}-x_{(1)}\log (1-\text{α})+x_{(1)}$

2.7 Bayes’ Theorem and Bayesian Statistics

The Law of Total Probability

If events A1,...,Ak partition a sample space S into mutually exclusive and exhaustive nonempty events, then the Law of Total Probability states that the total probability of an event B is given by

$\begin{array}{l}P\left(B\right)=P\left(A_{1}B\right)+P\left(A_{2}B\right)+...+P\left(A_{k}B\right)\\ =P\left(B|A_1\right)P\left(A_1\right)+P\left(B|A_2\right)P\left(A_2\right)+...+P\left(B|A_k\right)P\left(A_k\right)\\ ={\displaystyle \sum _{j=1}^{k}P\left(B|A_j\right)P\left(A_j\right).}\end{array}$

For continuous random variables X and Y we have the distributional form of the Law of Total Probability

$f\text{y}(y)={\displaystyle {\int }_{-\infty }^{\infty }f\text{y}|{\text{x}}_{\text{=}\infty }(y)f\text{x(}x\text{)}dx}$

For discrete random variables X and Y we can write the distributional form of the Law of Total Probability as

$f\text{y}(y)=P(Y=y){\displaystyle \sum _{x}P(Y_{=y}|{\text{X}}_{\text{=}x})}$

Bayes’ Theorem

Bayes’ Theorem provides a method for inverting conditional probabilities. In its simplest form, if A and B are events and P(B) > 0, then

$P(A|B)=\frac{P(B|A)P(A)}{P(B)}$

Often the Law of Total Probability is applied to compute P(B) in the denominator. These formulas follow from the definitions of conditional and joint probability.

For continuous random variables the distributional form of Bayes’ Theorem is

$f\text{x}|{\text{Y}}_{\text{=y}}(x)\frac{f\text{Y}|{\text{x}}_{\text{=}x}(y)f\text{x(}x\text{)}}{f\text{Y}(y)}=\frac{f\text{Y}|{\text{x}}_{\text{=}x}(y)f\text{x(}x\text{)}}{{\displaystyle {\int }_{-\infty }^{\infty }f\text{Y}|{\text{x}}_{\text{=}x}(y)f\text{x(}x\text{)}dx}}$

$f_{X|Y=y}(x)=P(X=x|\text{Y}=y)=\frac{P(Y=y|X=x)P(X=x)}{{\displaystyle {\sum }_x\left\{P(Y=y|X=x)P(X=x)\right\}}}$

These formulas follow from the definitions of conditional and joint probability.

Bayesian Statistics

In the frequentist approach to statistics, the parameters of a distribution are considered to be fixed but unknown constants. The Bayesian approach views the unknown parameters of a distribution as random variables. Thus, in Bayesian analysis, probabilities can be computed for parameters as well as the sample statistics.

Bayes’ Theorem allows one to revise his/her prior belief about an unknown parameter based on observed data. The prior belief reflects the relative weights that one assigns to the possible values for the parameters. Suppose that X has the density f(x|θ). The conditional density of θ given the sample observations x1, ... , xn is called the posterior density, defined by

$f_{\theta |x}(\theta )=\frac{f(x_1,\mathrm{...},x_n|\theta )f_{\theta }(\theta )}{{\displaystyle \int f(x_1,\mathrm{...},x_n|\theta )f_{\theta }(\theta )d\theta }}$

where fθ(θ) is the pdf of the prior distribution of θ. The posterior distribution summarizes our modified beliefs about the unknown parameters, taking into account the data that has been observed. Then one is interested in computing posterior quantities such as posterior means, posterior modes, posterior standard deviations, etc.

Note that any constant in the likelihood function cancels out of the posterior density. The basic relation is

$posterior\propto prior\times likelihood,$

which describes the shape of the posterior density up to a multiplicative constant. Often the evaluation of the constant is difficult and the integral cannot be obtained in closed form. However, Monte Carlo methods are available that do not require the evaluation of the constant in order to sample from the posterior distribution and estimate posterior quantities of interest. See e.g. [44, 103103, 106, 120, 228] on development of Markov Chain Monte Carlo sampling.

Readers are referred to Lee [171] for an introductory presentation of Bayesian statistics. Albert [5] is a good introduction to computational Bayesian methods with R. A textbook covering probability and mathematical statistics from both a classical and Bayesian perspective at an advanced undergraduate level is DeGroot and Schervish [64].

2.8 Markov Chains

In this section we briefly review discrete time, discrete state space Markov chains. A basic understanding of Markov chains is necessary background for Chapter 9 on Markov Chain Monte Carlo methods. Readers are referred to Ross [234, Ch. 4] for an excellent introduction to Markov chains.

A Markov chain is a stochastic process {Xt} indexed by time t ≥ 0. Our goal is to generate a chain by simulation, so we consider discrete time Markov chains. The time index will be the nonnegative integers, so that the process starts in state X0 and makes successive transitions to X1, X2,... , Xt, .... The set of possible values of Xt is the state space.

Suppose that the state space of a Markov chain is finite or countable. Without loss of generality, we can suppose that the states are 0, 1, 2, .... The sequence {Xt|t ≥ 0} is a Markov chain if

$P(X_{t+1}=j|X_0=i_0,X_1=i_1,\mathrm{...},X_{t-1}=i_{t-1},X_t=i)=P(X_{t+1}=j|X_t=i)$

for all pairs of states (i, j),t ≥ 0. In other words, the transition probability depends only on the current state, and not on the past.

If the state space is finite, the transition probabilitiesP(Xt+1|Xt) can be represented by a transition matrix] ℙ= (pij) where the entry pij is the probability that the chain makes a transition to state j in one step starting from state i. The probability that the chain moves from state i to state j in k steps is $p_{ij}^{(k)}$, and the Chapman-Kolmogorov equations (see e.g. [234, Ch. 4]) provide that the k-step transition probabilities are the entries of the matrix ℙk. That is, ℙ(k) = $p_{ij}^{(k)}$=ℙ(k) the kth power of the transition matrix.

A Markov chain is irreducible if all states communicate with all other states: given that the chain is in state i, there is a positive probability that the chain

can enter state j in finite time, for all pairs of states (i, j). A state i is recurrent if the chain returns to i with probability 1; otherwise state i is transient. If the expected time until the chain returns to i is finite, then i is nonnull or positive recurrent. The period of a state i is the greatest common divisor of the lengths of paths starting and ending at i. In an irreducible chain, the periods of all states are equal, and the chain is aperiodic if the states all have period 1. Positive recurrent, aperiodic states are ergodic. In a finite-state Markov chain all recurrent states are positive recurrent.

In an irreducible, ergodic Markov chain the transition probabilities converge to a stationary distribution π on the state space, independent of the initial state of the chain.

In a finite-state Markov chain, irreducibility and aperiodicity imply that for all states j

${\text{π}}_j=\underset{n\to \infty }{\lim }{p}_{ij}^{(n)}$

exists and is independent of the initial state i. The probability distribution π = {πj} is called the stationary distribution, and π is the unique nonnegative solution to the system of equations

$\begin{array}{cccc}{\text{π}}_j={\displaystyle \sum _{i=0}^{\infty }{\text{π}}_i{p}_{ij},}& j\ge 0;& {\displaystyle \sum _{j=0}^{\infty }{\pi }_j=1.}& (2.17)\end{array}$

We can interpret πj as the (limiting) proportion of time that the chain is in state j.

Example 2.7 (Finite state Markov chain)

Ross [234] gives the following example of a Markov chain model for mutations of DNA. A DNA nucleotide has four possible values. For each unit of time the model specifies that the nucleotide changes with probability 3α, for some0<α< 1/3. If it does change, then it is equally likely to change to any of the other three values. Thus pii = 1 − 3α and pij = 3α/3 = α, i≠j. If we number the states 1 to 4, the transition matrix is

$\begin{array}{cc}\mathbb{P}=\left[\begin{array}{llll}1-3\alpha \hfill & \alpha \hfill & \alpha \hfill & \alpha \hfill \\ \alpha \hfill & 1-3\alpha \hfill & \alpha \hfill & \alpha \hfill \\ \alpha \hfill & \alpha \hfill & 1-3\alpha \hfill & \alpha \hfill \\ \alpha \hfill & \alpha \hfill & \alpha \hfill & 1-3\alpha \hfill \end{array}\right]& (2.18)\end{array}$

where pij=ℙi,j is the probability of a mutation from state i to state j. The ith row of a transition matrix is the conditional probability distribution P(Xn+1 = j|Xn = i), j = 1, 2, 3, 4 of a transition to state j given that the process is currently in state i. Thus each row must sum to 1 (the matrix is row stochastic). This matrix happens to be doubly stochastic because the columns also sum to 1, but in general a transition matrix need only be row stochastic.

Suppose that α = 0.1. Then the two-step and the 16-step transition matrices are

$\begin{array}{cc}{\mathbb{P}}^2=\left[\begin{array}{llll}0.52\hfill & 0.16\hfill & 0.16\hfill & 0.16\hfill \\ 0.16\hfill & 0.52\hfill & 0.16\hfill & 0.16\hfill \\ 0.16\hfill & 0.16\hfill & 0.52\hfill & 0.16\hfill \\ 0.16\hfill & 0.16\hfill & 0.16\hfill & 0.52\hfill \end{array}\right]& {\mathbb{P}}^{16}=\left[\begin{array}{llll}0.2626\hfill & 0.2458\hfill & 0.2458\hfill & 0.2458\hfill \\ 0.2458\hfill & 0.2626\hfill & 0.2458\hfill & 0.2458\hfill \\ 0.2458\hfill & 0.2458\hfill & 0.2626\hfill & 0.2458\hfill \\ 0.2458\hfill & 0.2458\hfill & 0.2458\hfill & 0.2626\hfill \end{array}\right]\end{array}$

The three-step transition matrix is ℙ2ℙ = ℙ3, etc.The probability $p_{14}^{(2)}$ of transition from state 1 to state 4 in two steps is $p_{\mathrm{1,4}}^{(2)}=0.16$, and the probability that the process returns to state 2 from state 2 in 16 steps is $p_{22}^{(16)}={\mathbb{P}}_{\mathrm{2,2}}^{16}=0.2626$.

All entries of P are positive, hence all states communicate; the chain is irreducible and ergodic. The transition probabilities in every row are converging to the same stationary distribution π on the four states. The stationary distribution is the solution of equations (2.17); in this case $\pi (i)=\frac{1}{4}$,i = 1, 2, 3, 4. (In this example, it can be shown that the limiting probabilities do not depend on $\alpha :{\mathbb{P}}_{ii}^n=\frac{1}{4}+\frac{3}{4}{\left(1-4\alpha \right)}^n\to \frac{1}{4}asn\to \infty .)$

Example 2.8(Random walk)

An example of a discrete-time Markov chain with an infinite state space is the random walk. The state space is the set of all integers, and the transition probabilities are

$\begin{array}{lll}{p}_{i,i+1}=p,& & i=0\pm 1\pm 2,\mathrm{...},\\ p_{i,i-1}=1-p,& & i=0\pm 1\pm 2,\mathrm{...},\\ p_{i,j}& & j\notin \left\{i-1,i+1\right\}.\end{array}$

In the random walk model, at each transition a step of unit length is taken at random to the right with probability p or left with probability 1 − p. The state of the process at time n is the current location of the walker at time n. Another interpretation considers the gambler who bets $1 on a sequence of Bernoulli(p) trials and wins or loses $1 at each transition; if X0 = 0, the state of the process at time n is his gain or loss after n trials.

In the random walk model all states communicate, so the chain is irreducible. All states have period 2. For example, it is impossible to return to state 0 starting from 0 in an odd number of steps. The probability that the first return to 0 from state 0 occurs in exactly 2n steps is

$p_{00}^{(2n)}=\left(\begin{array}{c}2n\\ n\end{array}\right)p^n{\left(1-p\right)}^n=\frac{(2n)!}{n!n!}{\left(p\left(1-p\right)\right)}^n$

It can be shown that ${\displaystyle {\sum }_{n=1}^{\infty }{p}_{00}^{(2n)}<\infty }$ if and only if p≠ 1/2 Thus, the expected number of visits to 0 is finite if and only if p≠ 1/2. Recurrence and transience are class properties, hence the chain is recurrent if and only if p = 1/2 and otherwise all states are transient. When p = 1/2 the process is called a symmetric random walk. The symmetric random walk is discussed in Example 3.26.