Can anyone explain conjugate priors in simplest possible terms? I have been trying to understand the idea of conjugate priors in Bayesian statistics for a while but I simply don't get it. Can anyone explain the idea in the simplest possible terms, perhaps using the "Gaussian prior" as an example?
 A: I like using the notion of a "kernel" of a distribution. This is where you only leave in the parts that depend on the parameter. A few simple examples.
Normal kernel
$$p(\mu|a,b) = K^{-1} \times \exp(a\mu^2 +b\mu)$$
Where $K$ is the "normalising constant" $K=\int \exp(a\mu^2 +b\mu)d\mu=\sqrt{\frac{\pi}{-a}}\exp(-\frac{b^2}{4a})$
The connection with standard mean/variance parameters is $E(\mu|a,b)=-\frac{b}{2a}$ and $Var(\mu|a,b)=-\frac{1}{2a}$
Beta kernel
$$p(\theta|a,b)=K^{-1}\times \theta^a (1-\theta)^b$$
Where $K=\int \theta^a (1-\theta)^b d\theta = Beta(a+1,b+1)$
When we look at the likelihood function, we can do the same thing, and express it in "kernel form". For example with iid data
$$p(D|\mu)=\prod_{i=1}^n p(x_i|\mu)=Q\times f(\mu)$$
For some constant $Q$ and some function $f(\mu)$. If we can recognise this function as a kernel, then we can create a conjugate prior for that likelihood.
If we take the normal likelihood with unit variance, the above looks like
$$p(D|\mu)
=\prod_{i=1}^n p(x_i|\mu)
=\prod_{i=1}^n \frac{1}{\sqrt{2\pi}}\exp(-\frac{(x_i-\mu)^2}{2})
=\left[\prod_{i=1}^n \frac{1}{\sqrt{2\pi}}\right]\times \prod_{i=1}^n \exp(-\frac{(x_i-\mu)^2}{2})
=(2\pi)^{-\frac{n}{2}}\times\exp(-\sum_{i=1}^n\frac{(x_i-\mu)^2}{2})
=(2\pi)^{-\frac{n}{2}}\times\exp(-\sum_{i=1}^n\frac{x_i^2-2x_i\mu+\mu^2}{2})
=(2\pi)^{-\frac{n}{2}}\times\exp(-\sum_{i=1}^n\frac{x_i^2}{2})\times\exp(\mu\sum_{i=1}^n x_i-\mu^2\frac{n}{2})
=Q\times \exp(a\mu^2 +b\mu)$$
where $a=-\frac{n}{2}$ and $b=\sum_{i=1}^n x_i$ and $Q=(2\pi)^{-\frac{n}{2}}\times\exp(-\sum_{i=1}^n\frac{x_i^2}{2})$
This likelihood function has the same kernel as the normal distribution for $\mu$, so a conjugate prior for this likelihood is also the normal distribution.
$$p(\mu|a_0,b_0)=K_0^{-1}\exp(a_0\mu^2 +b_0\mu)$$
The posterior is then
$$p(\mu|D,a_0,b_0)\propto K_0^{-1}\exp(a_0\mu^2 +b_0\mu)\times Q\times \exp(a\mu^2 +b\mu)=K_0^{-1}\times Q\times \exp([a+a_0]\mu^2 +[b+b_0]\mu)\propto\exp([a+a_0]\mu^2 +[b+b_0]\mu)$$
Showing that the posterior is also a normal distribution, with updated parameters from the prior using the information in the data.
In some sense a conjugate prior acts similarly to adding "pseudo data" to the data observed, and then estimating the parameters.
A: A prior for a parameter will almost always have some specific functional form (written in terms of the density, generally). Let's say we restrict ourselves to one particular family of distributions, in which case choosing our prior reduces to choosing the parameters of that family.
For example, consider a normal model $Y_i \stackrel{_\text{iid}}{\sim} N(\mu,\sigma^2)$. For simplicity, let's also take $\sigma^2$ as known. This part of the model - the model for the data - determines the likelihood function.
To complete our Bayesian model, here we need a prior for $\mu$.
As mentioned above, commonly we might specify some distributional family for our prior for $\mu$ and then we only have to choose the parameters of that distribution (for example, often prior information may be fairly vague - like roughly where we want the probability to concentrate - rather than of very specific functional form, and we may have enough freedom to model what we want by choosing the parameters - say to match a prior mean and variance).
If it turns out that the posterior for $\mu$ is from the same family as the prior, then that prior is said to be "conjugate".
(What makes it turn out to be conjugate is the way it combines with the likelihood)
So in this case, let's take a Gaussian prior for $\mu$ (say $\mu\sim N(\theta,\tau^2)$). If we do that, we see that the posterior for $\mu$ is also Gaussian. Consequently, the Gaussian prior was a conjugate prior for our model above.
That's all there is to it really -- if the posterior is from the same family as the prior, it's a conjugate prior.
In simple cases you can identify a conjugate prior by inspection of the likelihood. For example, consider a binomial likelihood; dropping the constants, it looks like a beta density in $p$; and because of the way powers of $p$ and $(1-p)$ each combine, it will multiply by a beta prior to also give a product of powers of $p$ and $(1-p)$ ... so we can see immediately from the likelihood that the beta will be a conjugate prior for $p$ in the binomial likelihood.
In the Gaussian case it's easiest to see that it will happen by considering the log-densities and the log-likelihood; the log-likelihood will be quadratic in $\mu$ and the sum of two quadratics is quadratic, so a quadratic log-prior + quadratic log-likelihood gives a quadratic posterior (each of the coefficients of the highest order term will of course be negative).
A: For a given distribution family $D_{lik}$ of the likelihood (e.g. Bernoulli),
if the prior is of the same distribution family $D_{pri}$ as the posterior (e.g. Beta),
then $D_{pri}$ and $D_{lik}$ are conjugate distribution families and the prior is called a conjugate prior for the likelihood function.
Note:  $\underbrace{p(\theta|x)}_{\text{posterior}}
\sim 
\underbrace{p(x|\theta)}_{\text{likelihood}}
\cdot
\underbrace{p(\theta)}_{\text{prior}}$
A: If your model belongs to an exponential family, that is, if the density of the distribution is of the form
$$f(x|\theta)=h(x)\exp\{T(\theta)\cdot S(x)-\psi(\theta)\}\qquad
x\in\mathcal{X}\quad\theta\in\Theta$$
with respect to a given dominating measure (Lebesgue, counting, &tc.), where $t\cdot s$ denotes a scalar product over $\mathbb{R}^d$ and
$$T:\mathcal{X}\longrightarrow \mathbb{R}^d\qquad S:\Theta\longrightarrow \mathbb{R}^d$$
are measurable functions, the conjugate priors on $\theta$ are defined by densities of the form
$$\pi(\theta|\xi,\lambda)=C(\xi,\lambda)\exp\{T(\theta)\cdot \xi-\lambda\psi(\theta)\}$$
[with respect to an arbitrarily-chosen dominating measure $\text{d}\nu$ on $\Theta$] with
$$C(\xi,\lambda)^{-1}=\int_\Theta \exp\{T(\theta)\cdot \xi-\lambda\psi(\theta)\} \text{d}\nu<\infty$$
and $\lambda\in\Lambda\subset\mathbb{R}_+$, $\xi\in\Xi\subset \lambda T(\mathcal{X})$ 
The choice of the dominating measure is determinantal for the family of priors. If for instance one faces a Normal mean likelihood on $\mu$ as in Glen_b's answer, choosing the Lebesgue measure $\text{d}\mu$ as the dominating measure leads to Normal priors being conjugate. If instead one chooses $(1+\mu^2)^{-2}\text{d}\mu$ as the dominating measure, the conjugate priors are within the family of distributions with densities
$$\exp\{-\alpha(\mu-\mu_0)^2\} \qquad\alpha>0,\ \ \mu_0\in\mathbb R$$
with respect to this dominating measure, and are thus no longer Normal priors. This difficulty is essentially the same as the one of choosing a particular parameterisation of the likelihood and opting for the Lebesgue measure for this parameterisation. When faced with a likelihood function, there is no inherent (or intrinsic or reference) dominating measure on the parameter space.

Outside this exponential family setting, there is no non-trivial
  family of distributions with a fixed support that allows for conjugate
  priors. This is a consequence of the Darmois-Pitman-Koopman
  lemma.

