Conjugate prior explained

p(x \mid \theta) , the posterior distribution p(\theta \mid x) is in the same probability distribution family as the prior probability distribution p(\theta) , the prior and posterior are then called conjugate distributions with respect to that likelihood function and the prior is called a conjugate prior for the likelihood function p(x \mid \theta) .

A conjugate prior is an algebraic convenience, giving a closed-form expression for the posterior; otherwise, numerical integration may be necessary. Further, conjugate priors may give intuition by more transparently showing how a likelihood function updates a prior distribution.

The concept, as well as the term "conjugate prior", were introduced by Howard Raiffa and Robert Schlaifer in their work on Bayesian decision theory.^[1] A similar concept had been discovered independently by George Alfred Barnard.^[2]

Example

The form of the conjugate prior can generally be determined by inspection of the probability density or probability mass function of a distribution. For example, consider a random variable which consists of the number of successes

s in n Bernoulli trials with unknown probability of success q in [0,1]. This random variable will follow the binomial distribution, with a probability mass function of the formp(s) = {n \choose s}q ^{s} ( 1- q) ^{n-s}

The usual conjugate prior is the beta distribution with parameters (

\alpha , \beta ):p(q) = {q ^{\alpha-1} ( 1- q) ^{\beta-1} \over \Beta(\alpha,\beta)} where \alpha and \beta are chosen to reflect any existing belief or information (\alpha =1 and \beta =1 would give a uniform distribution) and \Beta(\alpha,\beta) is the Beta function acting as a normalising constant.

In this context,

\alpha and \beta are called hyperparameters (parameters of the prior), to distinguish them from parameters of the underlying model (here q ). A typical characteristic of conjugate priors is that the dimensionality of the hyperparameters is one greater than that of the parameters of the original distribution. If all parameters are scalar values, then there will be one more hyperparameter than parameter; but this also applies to vector-valued and matrix-valued parameters. (See the general article on the exponential family, and also consider the Wishart distribution, conjugate prior of the covariance matrix of a multivariate normal distribution, for an example where a large dimensionality is involved.)

If we sample this random variable and get

s successes and f = n - s failures, then we have\begin{align} P(s, f \mid q = x) & = {s + f \choose s} x ^{s(1-x)} ^{f,\\} P(q = x) & = {x ^{\alpha-1} ( 1- x) ^{\beta-1} \over \Beta(\alpha,\beta)},\\ P(q = x \mid s,f) & = (P(s, f \mid x)P(x))/(\int P(s, f \mid y)P(y)dy) \\ & = {{{s + f \choose s} x ^{s+\alpha-1} ( 1- x) ^{f+\beta-1} / \Beta(\alpha,\beta)} \over \int ^{1} _{y =0} \left({s + f \choose s} y ^{s+\alpha-1} ( 1- y) ^{f+\beta-1} / \Beta(\alpha,\beta)\right) dy} \\ & = {x ^{s+\alpha-1} ( 1- x) ^{f+\beta-1} \over \Beta(s + \alpha,f + \beta)}, \end{align}

which is another Beta distribution with parameters

(\alpha + s, \beta + f) . This posterior distribution could then be used as the prior for more samples, with the hyperparameters simply adding each extra piece of information as it comes.

Interpretations

Pseudo-observations

It is often useful to think of the hyperparameters of a conjugate prior distribution corresponding to having observed a certain number of pseudo-observations with properties specified by the parameters. For example, the values

\alpha and \beta of a beta distribution can be thought of as corresponding to \alpha -1 successes and \beta -1 failures if the posterior mode is used to choose an optimal parameter setting, or \alpha successes and \beta failures if the posterior mean is used to choose an optimal parameter setting. In general, for nearly all conjugate prior distributions, the hyperparameters can be interpreted in terms of pseudo-observations. This can help provide intuition behind the often messy update equations and help choose reasonable hyperparameters for a prior.

Dynamical system

One can think of conditioning on conjugate priors as defining a kind of (discrete time) dynamical system: from a given set of hyperparameters, incoming data updates these hyperparameters, so one can see the change in hyperparameters as a kind of "time evolution" of the system, corresponding to "learning". Starting at different points yields different flows over time. This is again analogous with the dynamical system defined by a linear operator, but note that since different samples lead to different inferences, this is not simply dependent on time but rather on data over time. For related approaches, see Recursive Bayesian estimation and Data assimilation.

Practical example

Suppose a rental car service operates in your city. Drivers can drop off and pick up cars anywhere inside the city limits. You can find and rent cars using an app.

Suppose you wish to find the probability that you can find a rental car within a short distance of your home address at any time of day.

Over three days you look at the app and find the following number of cars within a short distance of your home address:

x = [ 3 , 4 , 1 ]

Suppose we assume the data comes from a Poisson distribution. In that case, we can compute the maximum likelihood estimate of the parameters of the model, which is $\backslash lambda\; =\; \backslash frac\; \backslash approx\; 2.67.$ Using this maximum likelihood estimate, we can compute the probability that there will be at least one car available on a given day: $p(x>0\; |\; \backslash lambda\; \backslash approx\; 2.67)\; =\; 1\; -\; p(x=0\; |\; \backslash lambda\; \backslash approx\; 2.67)\; =\; 1-\backslash frac\; \backslash approx\; 0.93$

This is the Poisson distribution that is the most likely to have generated the observed data

x . But the data could also have come from another Poisson distribution, e.g., one with λ = 3 , or λ = 2 , etc. In fact, there is an infinite number of Poisson distributions that could have generated the observed data. With relatively few data points, we should be quite uncertain about which exact Poisson distribution generated this data. Intuitively we should instead take a weighted average of the probability of p(x> 0 | λ ) for each of those Poisson distributions, weighted by how likely they each are, given the data we've observed x .p(x| x ) = \int _{\theta} p(x|\theta)p(\theta| x )d\theta, where x is a new data point, x is the observed data and \theta are the parameters of the model. Using Bayes' theorem we can expand p(\theta| x ) = (p( x |\theta)p(\theta))/(p( x )) , therefore p(x| x ) = \int _{\theta} p(x|\theta) (p( x |\theta)p(\theta))/(p( x )) d\theta. Generally, this integral is hard to compute. However, if you choose a conjugate prior distribution p(\theta) , a closed-form expression can be derived. This is the posterior predictive column in the tables below.

Returning to our example, if we pick the Gamma distribution as our prior distribution over the rate of the Poisson distributions, then the posterior predictive is the negative binomial distribution, as can be seen from the table below. The Gamma distribution is parameterized by two hyperparameters

\alpha, \beta , which we have to choose. By looking at plots of the gamma distribution, we pick \alpha = \beta = 2 , which seems to be a reasonable prior for the average number of cars. The choice of prior hyperparameters is inherently subjective and based on prior knowledge.

Given the prior hyperparameters

\alpha and \beta we can compute the posterior hyperparameters $\backslash alpha\text{'}\; =\; \backslash alpha\; +\; \backslash sum\_i\; x\_i\; =\; 2\; +\; 3+4+1\; =\; 10$ and $\backslash beta\text{'}\; =\; \backslash beta\; +\; n\; =\; 2+3\; =\; 5$

Given the posterior hyperparameters, we can finally compute the posterior predictive of $p(x>0|\backslash mathbf)\; =\; 1-p(x=0|\backslash mathbf)\; =\; 1\; -\; NB\backslash left(0\backslash ,\; |\backslash ,\; 10,\; \backslash frac\backslash right)\; \backslash approx\; 0.84$

This much more conservative estimate reflects the uncertainty in the model parameters, which the posterior predictive takes into account.

Table of conjugate distributions

Let n denote the number of observations. In all cases below, the data is assumed to consist of n points

x _{1,\ldots,x} _{n} (which will be random vectors in the multivariate cases).

If the likelihood function belongs to the exponential family, then a conjugate prior exists, often also in the exponential family; see Exponential family: Conjugate distributions.

When the likelihood function is a discrete distribution

Likelihood p(x _{i}	\theta)!	Model parameters \theta	Conjugate prior (and posterior) distribution p(\theta	\Theta), p(\theta	\mathbf,\Theta) = p(\theta	\Theta') !	Prior hyperparameters \Theta	Posterior hyperparameters \Theta'	Interpretation of hyperparameters	Posterior predictive p(\tilde{x}	\mathbf, \Theta) = p(\tilde	\Theta')
	p (probability)		\alpha, \beta\in R	\alpha + \sum ^{n} _{i =1} x _{i,} \beta + n - \sum ^{n} _{i =1} x _{i}	\alpha successes, \beta failures	p(\tilde{x} =1 ) = (\alpha')/(\alpha' + \beta') (Bernoulli)
Binomial with known number of trials, m	p (probability)		\alpha, \beta\in R	\alpha + \sum ^{n} _{i =1} x _{i,} \beta + \sum ^{nN} _{i} - \sum ^{n} _{i =1} x _{i}	\alpha successes, \beta failures	\operatorname{BetaBin}(\tilde{x}	\alpha',\beta') (beta-binomial)
Negative binomial with known failure number, r	p (probability)		\alpha, \beta\in R	\alpha + rn , \beta + \sum ^{n} _{i =1} x _{i}	\alpha total successes, \beta failures (i.e., (\beta)/(r) experiments, assuming r stays fixed)	\operatorname{BetaNegBin}(\tilde{x}	\alpha',\beta')(beta-negative binomial)
Poisson	λ (rate)	Gamma	k, \theta\in R	k + \sum ^{n} _{i =1} x _{i, (\theta)/(n \theta +1)}	k total occurrences in (1)/(\theta) intervals	\operatorname{NB}\left(\tilde{x}\mid k', (1)/(\theta' +1) \right) (negative binomial)
			\alpha, \beta	\alpha + \sum ^{n} _{i =1} x _{i} , \beta + n	\alpha total occurrences in \beta intervals	\operatorname{NB}\left(\tilde{x}\mid\alpha', (\beta')/(1+ \beta') \right) (negative binomial)
	p (probability vector), k (number of categories; i.e., size of p)		\boldsymbol\alpha\in R ^{k}	\boldsymbol\alpha + (c _{1,} \ldots, c _{k),} where c _{i} is the number of observations in category i	\alpha _{i} occurrences of category i	\begin{align} p(\tilde{x} = i) & = ({\alpha _{i})/('}{\sum) _{i} {\alpha _{i}'}} \\ & = (\alpha _{i} + c _{i})/(\sum _{i} \alpha _{i} + n) \end{align} (categorical)
	p (probability vector), k (number of categories; i.e., size of p)		\boldsymbol\alpha\in R ^{k}	\boldsymbol\alpha + \sum ^{n x} _{i}	\alpha _{i} occurrences of category i	\operatorname{DirMult}(\tilde{ x }\mid\boldsymbol\alpha') (Dirichlet-multinomial)
Hypergeometric with known total population size, N	M (number of target members)		n = N, \alpha, \beta	\alpha + \sum ^{n} _{i =1} x _{i,} \beta + \sum ^{nN} _{i} - \sum ^{n} _{i =1} x _{i}	\alpha successes, \beta failures
	p0 (probability)		\alpha, \beta\in R	\alpha + n, \beta + \sum ^{n} _{i =1} x _{i}	\alpha experiments, \beta total failures

When likelihood function is a continuous distribution

Likelihood p(x _{i}	\theta)!	Model parameters \theta	Conjugate prior (and posterior) distribution p(\theta	\Theta), p(\theta	\mathbf,\Theta) = p(\theta	\Theta') !	Prior hyperparameters \Theta	Posterior hyperparameters \Theta'	Interpretation of hyperparameters	Posterior predictive p(\tilde{x}	\mathbf, \Theta) = p(\tilde	\Theta')
Normal with known variance σ^2	μ (mean)		\mu _{0,} \sigma ^{2} _{0}	(1)/(\left( (1)/(+ (n)/(\sigma ^{2}) \sigma ^{2} _{0})) (\mu _{0})/(\sigma ^{2} _{0}) + (\sum ^{n} _{i =1} x _{i})/(\right), \left( \sigma ^{2}) (1)/(\sigma ^{2} _{0}) + (n)/(\sigma ^{2}) \right) ^{-1}	mean was estimated from observations with total precision (sum of all individual precisions) 1 /\sigma ^{2} _{0} and with sample mean \mu _{0}	l{N}(\tilde{x}	\mu_0', ' +\sigma^2)
Normal with known precision τ	μ (mean)		\mu _{0,} \tau ^{-1} _{0}	(\tau _{0} \mu _{0} + \tau \sum ^{n} _{i =1} x _{i})/(\tau _{0} + n \tau) , \left(\tau _{0} + n \tau\right) ^{-1}	mean was estimated from observations with total precision (sum of all individual precisions)\tau _{0} and with sample mean \mu _{0}	l{N}\left(\tilde{x}\mid\mu _{0',} (1)/(+ \tau _{0'}) (1)/(\tau) \right)
Normal with known mean μ	σ^2 (variance)		\alpha,\beta	\alpha + (n)/(2) , \beta + (\sum ^{2} _{i - \mu)})/() }	variance was estimated from 2 \alpha observations with sample variance \beta/\alpha (i.e. with sum of squared deviations 2 \beta , where deviations are from known mean \mu )	t _{2\alpha'} (\tilde{x}	\mu,\sigma^2 = \beta'/\alpha')
Normal with known mean μ	σ^2 (variance)		\nu, \sigma ^{2} _{0}	\nu + n, (\nu\sigma ^{2} _{0} + \sum ^{n} _{i =1} (x ^{2} _{i - \mu)})/(\nu + n)	variance was estimated from \nu observations with sample variance \sigma ^{2} _{0}	t _{\nu'} (\tilde{x}	\mu,')
Normal with known mean μ	τ (precision)		\alpha, \beta	\alpha + (n)/(2) , \beta + (\sum ^{n} _{i =1} (x ^{2} _{i - \mu)})/(2)	precision was estimated from 2 \alpha observations with sample variance \beta/\alpha (i.e. with sum of squared deviations 2 \beta , where deviations are from known mean \mu )	t _{2\alpha'} (\tilde{x}\mid\mu,\sigma ^{2} = \beta'/\alpha')
Normal^[3]	μ and σ^2 Assuming exchangeability	Normal-inverse gamma	\mu _{0} , \nu , \alpha , \beta	(\nu\mu _{0+n\bar{x})/() },\, \nu+n,\, \alpha+\frac,\, \beta + \tfrac{ 1 }{ 2 } \sum ^{n} _{i =1} (x _{i} - \bar{x}) ^{2} + (n\nu)/(\nu + n) ((\bar{x)/(- \mu) ^{2 }{ 2 }} ^{0 )} \bar{x} is the sample mean	mean was estimated from \nu observations with sample mean \mu _{0} ; variance was estimated from 2 \alpha observations with sample mean \mu _{0} and sum of squared deviations 2 \beta	t _{2\alpha'} \left(\tilde{x}\mid\mu', (\beta'(\nu' +1 ))/(\nu' \alpha') \right)
	μ and τ Assuming exchangeability	Normal-gamma	\mu _{0} , \nu , \alpha , \beta	(\nu\mu _{0+n\bar{x})/() },\, \nu+n,\, \alpha+\frac,\, \beta + \tfrac{ 1 }{ 2 } \sum ^{n} _{i =1} (x _{i} - \bar{x}) ^{2} + (n\nu)/(\nu + n) ((\bar{x)/(- \mu) ^{2 }{ 2 }} ^{0 )} \bar{x} is the sample mean	mean was estimated from \nu observations with sample mean \mu _{0} , and precision was estimated from 2 \alpha observations with sample mean \mu _{0} and sum of squared deviations 2 \beta	t _{2\alpha'} \left(\tilde{x}\mid\mu', (\beta'(\nu' +1 ))/(\alpha'\nu') \right)
Multivariate normal with known covariance matrix Σ	μ (mean vector)		\boldsymbol{\boldsymbol\mu} _{0,} \boldsymbol\Sigma _{0}	\left(\boldsymbol\Sigma ^{-1} _{0} + n\boldsymbol\Sigma ^{-1} \right) ^{-1} \left( \boldsymbol\Sigma ^{-1} _{0} \boldsymbol\mu _{0} + n \boldsymbol\Sigma ^{-1} \bar{x } \right), \left(\boldsymbol\Sigma ^{-1} _{0} + n\boldsymbol\Sigma ^{-1} \right) ^{-1} \bar{x } is the sample mean	mean was estimated from observations with total precision (sum of all individual precisions)\boldsymbol\Sigma ^{-1} _{0} and with sample mean \boldsymbol\mu _{0}	l{N}(\tilde{ x }\mid', ' +\boldsymbol\Sigma)
Multivariate normal with known precision matrix Λ	μ (mean vector)		\boldsymbol\mu _{0,} \boldsymbol Λ _{0}	\left(\boldsymbol Λ _{0} + n\boldsymbol Λ \right) ^{-1} \left( \boldsymbol Λ _{0\boldsymbol\mu} _{0} + n \boldsymbol Λ \bar{x } \right),\, \left(\boldsymbol\Lambda_0 + n\boldsymbol\Lambda\right)\bar{x } is the sample mean	mean was estimated from observations with total precision (sum of all individual precisions)\boldsymbol Λ _{0} and with sample mean \boldsymbol\mu _{0}	l{N}\left(\tilde{ x }\mid', ^ + \boldsymbol\Lambda^\right)
Multivariate normal with known mean μ	Σ (covariance matrix)		\nu , \boldsymbol\Psi	n + \nu , \boldsymbol\Psi + \sum ^{n} _{i =1} ( x _{i} - \boldsymbol\mu) ( x _{i} - \boldsymbol\mu) ^{T}	covariance matrix was estimated from \nu observations with sum of pairwise deviation products \boldsymbol\Psi	t _{\nu'-p+1} \left(\tilde{ x }	\boldsymbol\mu,\frac\boldsymbol\Psi'\right)
Multivariate normal with known mean μ	Λ (precision matrix)		\nu , V	n + \nu , \left( V ^{-1} + \sum ^{n} _{i =1} ( x _{i} - \boldsymbol\mu) ( x _{i} - \boldsymbol\mu) ^{T\right)} ^{-1}	covariance matrix was estimated from \nu observations with sum of pairwise deviation products V ^{-1}	t _{\nu'-p+1} \left(\tilde{ x }\mid\boldsymbol\mu,\frac^\right)
	μ (mean vector) and Σ (covariance matrix)		\boldsymbol\mu _{0} , \kappa _{0} , \nu _{0} , \boldsymbol\Psi	(\kappa _{0\boldsymbol\mu} _{0+n\bar{x})/() },\, \kappa_0+n,\, \nu_0+n,\, \boldsymbol\Psi + C + (\kappa _{0} n)/(\kappa _{0+n}) ( \bar{x }-\boldsymbol\mu_0)(\mathbf-\boldsymbol\mu_0)^T \bar{x } is the sample meanC = \sum ^{n} _{i =1} ( x _{i} - \bar{x }) (\mathbf - \mathbf)^T	mean was estimated from \kappa _{0} observations with sample mean \boldsymbol\mu _{0} ; covariance matrix was estimated from \nu _{0} observations with sample mean \boldsymbol\mu _{0} and with sum of pairwise deviation products \boldsymbol\Psi = \nu _{0\boldsymbol\Sigma} _{0}	t _{{\nu _{0}} ' - p +1 }\left(\tilde{ x }	',\frac\boldsymbol\Psi'\right)
	μ (mean vector) and Λ (precision matrix)		\boldsymbol\mu _{0} , \kappa _{0} , \nu _{0} , V	(\kappa _{0\boldsymbol\mu} _{0+n\bar{x})/() },\, \kappa_0+n,\, \nu_0+n,\, \left( V ^{-1} + C + (\kappa _{0} n)/(\kappa _{0+n}) ( \bar{x }-\boldsymbol\mu_0)(\mathbf-\boldsymbol\mu_0)^T\right)^ \bar{x } is the sample meanC = \sum ^{n} _{i =1} ( x _{i} - \bar{x }) (\mathbf - \mathbf)^T	mean was estimated from \kappa _{0} observations with sample mean \boldsymbol\mu _{0} ; covariance matrix was estimated from \nu _{0} observations with sample mean \boldsymbol\mu _{0} and with sum of pairwise deviation products V ^{-1}	t _{{\nu _{0}} ' - p +1 }\left(\tilde{ x }\mid ', \frac^\right)
	U( 0 ,\theta)		x _{m} , k	max \{x _{1,\ldots,x} _{n,x} _{m\},} k + n	k observations with maximum value x _{m}
Pareto with known minimum xm	k (shape)		\alpha, \beta	\alpha + n, \beta + \sum ^{n} _{i =1} ln (x _{i})/(x _{m})	\alpha observations with sum \beta of the order of magnitude of each observation (i.e. the logarithm of the ratio of each observation to the minimum x _{m} )
Weibull with known shape β	θ (scale)	Inverse gamma	a, b	a + n, b + \sum ^{n} _{i =1} x ^{\beta} _{i}	a observations with sum b of the βth power of each observation
Log-normal	Same as for the normal distribution after applying the natural logarithm to the data for the posterior hyperparameters. Please refer to to see the details.
	λ (rate)		\alpha, \beta	\alpha + n, \beta + \sum ^{n} _{i =1} x _{i}	\alpha observations that sum to \beta ^[4]	\operatorname{Lomax}(\tilde{x}\mid\beta',\alpha') (Lomax distribution)
Gamma with known shape α	β (rate)		\alpha _{0,} \beta _{0}	\alpha _{0+n\alpha,} \beta _{0+\sum} ^{n} ^{i =1} x _{i}	\alpha _{0/\alpha} observations with sum \beta _{0}	\operatorname{CG}(\tilde{ x }\mid\alpha,',')=\operatorname(\tilde	\alpha,',1,')
Inverse Gamma with known shape α	β (inverse scale)		\alpha _{0,} \beta _{0}	\alpha _{0+n\alpha,} \beta _{0+\sum} ^{n} ^{i =1} (1)/(x _{i})	\alpha _{0/\alpha} observations with sum \beta _{0}
Gamma with known rate β	α (shape)	\propto (a ^{\alpha-1} \beta ^{\alpha})/(\Gamma(\alpha) ^{b})	a, b, c	a \prod ^{n} _{i =1} x _{i,} b + n, c + n	b or c observations (b for estimating \alpha , c for estimating \beta ) with product a
Gamma	α (shape), β (inverse scale)	\propto (p ^{\alpha-1} e ^{-\beta})/(\Gamma(\alpha) ^{r} \beta ^{-\alpha})	p, q, r, s	p \prod ^{n} _{i =1} x _{i,} q + \sum ^{n} _{i =1} x _{i,} r + n, s + n	\alpha was estimated from r observations with product p ; \beta was estimated from s observations with sum q
Beta	α, β	\propto (\Gamma(\alpha + \beta) ^{k} p ^{\alpha} q ^{\beta})/(\Gamma(\alpha) ^{k\Gamma(\beta)} ^{k})	p, q, k	p \prod ^{n} _{i =1} x _{i,} q \prod ^{n} _{i =1} ( 1- x _{i),} k + n	\alpha and \beta were estimated from k observations with product p and product of the complements q

See also

• Beta-binomial distribution

Notes and References

1. [Howard Raiffa]
2. Jeff Miller et al. Earliest Known Uses of Some of the Words of Mathematics, "conjugate prior distributions". Electronic document, revision of November 13, 2005, retrieved December 2, 2005.
3. A different conjugate prior for unknown mean and variance, but with a fixed, linear relationship between them, is found in the normal variance-mean mixture, with the generalized inverse Gaussian as conjugate mixing distribution.
4. Book: Liu . Han . Statistical Machine Learning . Wasserman . Larry . 2014 . 314.