of Basics Bayesian Inference
A frequentist thinks of unknown parameters as fixed A Bayesian thinks of parameters as random, and thus having distributions (just like the data). A Bayesian writes down a prior guess for , and combines it with the likelihood for the observed data Y to obtain the posterior distribution of . All statistical inferences then follow from summarizing the posterior. This approach expands the class of candidate models, and facilitates hierarchical modeling, where it is important to properly account for various sources of uncertainty (e.g. spatial vs. nonspatial heterogeneity) The classical (frequentist) approach to estimation is not "wrong", but it is "limited in scope"!
Chapter 2: Basics of Point-Referenced Data Models p. 1/2
Basics of Bayesian Inference
As usual, we start with a model f (y|) for the observed data y = (y1 , . . . , yn ) given the unknown parameters = (1 , . . . , K ) Add a prior distribution (|), where is a vector of hyperparameters. The posterior distribution for is given by
p(y, |) = p(|y, ) = p(y|) = p(y, |) p(y, |) d
f (y|)(|) f (y|)(|) . = m(y|) f (y|)(|) d
We refer to this formula as Bayes' Theorem.
Chapter 2: Basics of Point-Referenced Data Models p. 2/2
Basics of Bayesian Inference
When is not known, a second stage (hyperprior) distribution h() will often be required, so that
p(y, ) = p(|y) = p(y) f (y|)(|)h() d . f (y|)(|)h() dd
Alternatively, we might replace in p(|y, ) by an ^ estimate ; this is called empirical Bayes analysis Note that posterior information prior information 0 , with the second "" replaced by "=" only if the prior is noninformative (which is often uniform, or "flat").
Chapter 2: Basics of Point-Referenced Data Models p. 3/2
Illustration of Bayes' Theorem
Suppose f (y|) = N (y|, 2 ), and > 0 known If we take (|) = N (|, 2 ) where = (, ) is fixed and known, then it is easy to show that
p(|y) = N 2 2 2 2 2 + 2 y, 2 2 2 + + + 2 .
Note that The posterior mean E(|y) is a weighted average of the prior mean and the data value y , with weights depending on our relative uncertainty the posterior precision (reciprocal of the variance) is equal to 1/ 2 + 1/ 2 , which is the sum of the likelihood and prior precisions.
Chapter 2: Basics of Point-Referenced Data Models p. 4/2
Illustration (continued)
As a concrete example, let = 2, = 1, y = 6, and = 1:
prior posterior with n = 1 posterior with n = 10
1.2 0.0 0.2 0.4 0.6 0.8 1.0
-2
0
2
4
6
8
When n = 1, prior and likelihood receive equal weight When n = 10, the data dominate the prior The posterior variance goes to zero as n
Chapter 2: Basics of Point-Referenced Data Models p. 5/2
Notes on priors
The prior here is conjugate: it leads to a posterior distribution for that is available in closed form, and is a member of the same distributional family as the prior. Note that setting 2 = corresponds to an arbitrarily vague (or noninformative) prior. The posterior is then
p (|y) = N |y, 2 /n ,
the same as the likelihood! The limit of the conjugate (normal) prior here is a uniform (or "flat") prior, and thus the posterior is the renormalized likelihood. The flat prior is appealing but improper here, since p()d = +. However, the posterior is still well defined, and so improper priors are often used!
Chapter 2: Basics of Point-Referenced Data Models p. 6/2
A more general example
Let Y be an n 1 data vector, X an n p matrix of covariates, and adopt the likelihood and prior structure,
Y| Nn (X, ) and Np (A, V )
Then the posterior distribution of |Y is
|Y N (Dd, D) , where D-1 = X T -1 X + V -1 and d = X T -1 Y + V -1 A. V -1 = 0 delivers a "flat" prior; if = 2 Ip , we get ^ |Y N , 2 (X X)-1 , where
^ = (X X)-1 X y usual likelihood approach!
Chapter 2: Basics of Point-Referenced Data Models p. 7/2
Bayesian estimation
Choose an appropriate measure of centrality: the posterior mean, median, or mode.
Point estimation:
Consider qL and qU , the /2- and (1 - /2)-quantiles of p(|y):
Interval estimation:
qL
p(|y)d = /2 and
- qU
p(|y)d = 1 - /2 .
Then clearly P (qL < < qU |y) = 1 - ; our confidence that lies in (qL , qU ) is 100 (1 - )%. Thus this interval is a 100 (1 - )% credible set ("Bayesian CI") for . This interval is relatively easy to compute, and enjoys a direct interpretation ("The probability that lies in (qL , qU ) is (1 - )") that the frequentist interval does not.
Chapter 2: Basics of Point-Referenced Data Models p. 8/2
Ex: Y Bin(10, ), U (0, 1), yobs = 7
2.0 0.0
0.0
0.5
1.0
1.5
2.5
3.0
0.2
0.4
0.6
0.8
1.0
Plot Beta(yobs + 1, n - yobs + 1) = Beta(8, 4) posterior in R/S: > theta <- seq(from=0, to=1, length=101) > yobs <- 7; n <- 10 > plot(theta, dbeta(theta, yobs+1, n-yobs+1), type="l") Add 95% equal-tail Bayesian CI (dotted vertical lines): > abline(v=qbeta(.5, yobs+1, n-yobs+1)) > abline(v=qbeta(c(.025, .975), yobs+1, n-yobs+1), lty=2)
Chapter 2: Basics of Point-Referenced Data Models p. 9/2
Bayesian hypothesis testing
Classical approach bases accept/reject decision on p-value = P {T (Y) more "extreme" than T (yobs )|, H0 } , where "extremeness" is in the direction of HA Several troubles with this approach: hypotheses must be nested p-value can only offer evidence against the null p-value is not the "probability that H0 is true" (but is often erroneously interpreted this way) As a result of the dependence on "more extreme" T (Y) values, two experiments with different designs but identical likelihoods could result in different p-values, violating the Likelihood Principle
Chapter 2: Basics of Point-Referenced Data Models p. 10/2
Bayesian hypothesis testing (cont'd)
Bayesian approach: Select the model with the largest posterior probability, P (Mi |y) = p(y|Mi )p(Mi )/p(y), where p(y|Mi ) =
f (y| i , Mi )i ( i )d i .
For two models, the quantity commonly used to summarize these results is the Bayes factor,
P (M1 |y)/P (M2 |y) p(y | M1 ) BF = = , P (M1 )/P (M2 ) p(y | M2 )
i.e., the likelihood ratio if both hypotheses are simple Problem: If i ( i ) is improper, then p(y|Mi ) necessarily is as well = BF is not well-defined!...
Chapter 2: Basics of Point-Referenced Data Models p. 11/2
Bayesian testing hypothesis via DIC
A generalization of the Akaike Information Criterion (AIC) to the case of hierarchical models based on the posterior distribution of the deviance statistic,
D() = -2 log f (y|) + 2 log h(y) ,
where f (y|) is the likelihood and h(y) is any standardizing function of the data alone Summarize the fit of a model by the posterior expectation of the deviance, D = E|y [D] Summarize the complexity of a model by the effective number of parameters,
pD = E|y [D] - D(E|y []) = D - D() .
Chapter 2: Basics of Point-Referenced Data Models p. 12/2
Bayesian hypothesis testing via DIC
The Deviance Information Criterion (DIC) is then
DIC = D + pD = 2D - D() ,
with smaller values indicating preferred models. Both building blocks of DIC and pD , E|y [D] and D(E|y []), are easily estimated via MCMC methods, and in fact are automatic within WinBUGS. While pD has a scale (effective model size), DIC does not, so only differences in DIC across models matter. DIC can be sensitive to parametrization and "focus" (i.e., what is considered to be part of the likelihood) f (y|): "focused on " p(y|) = f (y|)p(|)d : "focused on "
Chapter 2: Basics of Point-Referenced Data Models p. 13/2
Bayesian computation
Recall Monte Carlo integration: Suppose h(), and we seek E[f ()] =
1 = ^ N
N j=1 w.p. 1
f ()h()d. Then if i h(),
iid
f (j ) E[f ()] as N ( SLLN )
Since is itself a sample mean of independent ^ observations, V ar(^ ) = V ar[f ()]/N . This means
se(^ ) = ^ 1 N (N - 1)
N
[f (j ) - ]2 . ^
j=1
= use to get 95% (frequentist!) CI for
Chapter 2: Basics of Point-Referenced Data Models p. 14/2
Gibbs sampling
Suppose the joint distribution of = (1 , . . . , K ) is uniquely determined by the full conditional distributions, {pi (i |j=i ), i = 1, . . . , K}.
Gibbs Sampler:
Given an arbitrary set of starting values Draw 1
(1) (1) (0)
(0) (0) {1 , . . . , K },
p1 (1 |2 , . . . , K ), p2 (2 |1 , 3 , . . . , K ), . . . pK (K |1 , . . . , K-1 ),
(1) (1) (1) (0) (0)
(0)
Draw 2
Draw K
(1)
Under mild conditions,
(1 , . . . , K ) (1 , , K ) p as t .
(t) (t) d
Chapter 2: Basics of Point-Referenced Data Models p. 15/2
Gibbs sampling (cont'd)
For t sufficiently large (say, bigger than t0 ), { (t) }T 0 +1 t=t is a (correlated) sample from the true posterior. We might therefore use a sample mean to estimate the posterior mean, i.e.,
1 E(i |y) = T - t0
T
i .
t=t0 +1
(t)
The time from t = 0 to t = t0 is commonly known as the burn-in period; one can safely adapt (change) an MCMC algorithm during this preconvergence period, since these samples will be discarded anyway
Chapter 2: Basics of Point-Referenced Data Models p. 16/2
Gibbs sampling (cont'd)
In practice, we may actually run m parallel Gibbs sampling chains, instead of only 1, for some modest m (say, m = 5). Discarding the burn-in period, we obtain
1 E(i |y) = m(T - t0 )
m T (t) i,j j=1 t=t0 +1
,
where now the j subscript indicates chain number. What happens if the full conditional p(i |j=i , y) is not available in closed form? Typically, p(i |j=i , y) will be available up to proportionality constant, since it is proportional to the portion of f (y|)p() that involves i .
Chapter 2: Basics of Point-Referenced Data Models p. 17/2
Metropolis algorithm
Suppose the true joint posterior for has density p(u) wrt some measure . Choose an candidate density q( | (t-1) ) that is a valid density function for every possible value of the conditioning variable (t-1) , and satisfies q( | (t-1) ) = q( (t-1) | ), i.e., q is symmetric in its arguments. Given a starting value (0) at iteration t = 0, the algorithm proceeds as follows:
Chapter 2: Basics of Point-Referenced Data Models p. 18/2
Metropolis algorithm (cont'd)
Metropolis Algorithm:
For (t 1 ...