Maximum likelihood and Bayesian estimators for Bayesian networks

Motivation

A Bayesian network is a graphical model that captures the conditional dependence of variables via a directed acyclic graph. Consider, for example, the following mixture model in which the outcome of a (possibly unfair) coin toss determines which of two normal distributions to sample from.

Latent (a.k.a. unobservable) and observable variables are distinguished in the graphical model by the use of rectangular and elliptic nodes. Only the final result, $Y$, is observable. The coin toss, for example, is not. How does one estimate the parameters $\pi, \mu_0, \sigma_0, \mu_1, \sigma_1$ of this model?

Using the mixture model as a running example, this expository post considers maximum likelihood estimators and Bayesian estimators (e.g., maximum a posteriori estimators and credible intervals) to answer the above. More often than not, Bayesian networks do not have closed forms for these estimators (or they are unwieldy) and hence iterative algorithms are needed. Specifically, this post introduces classical optimization methods and expectation maximization (EM) for computing an MLE and Markov Chain Monte Carlo (MCMC) for computing Bayesian estimators.

Transformation to unconstrained space

The first thing to note is that some of the parameters are constrained. Namely, $\pi \in (0, 1)$ and $\sigma_i \in (0, \infty)$. Some of the methods described in this post are either incompatible with or perform poorly in the presence of constraints. In light of this, given parameters $\pi, \mu_0, \sigma_0, \mu_1, \sigma_1$, let

Note, in particular, that $\theta$ takes values in all of $\mathbb{R}^d$ with $d=5$.

def unconstrain(π, μ_0, σ_0, μ_1, σ_1):
    """Transform to unconstrained space."""
    logit_π = np.log(π / (1. - π))
    log_σ_0 = np.log(σ_0)
    log_σ_1 = np.log(σ_1)
    return logit_π, μ_0, log_σ_0, μ_1, log_σ_1

def constrain(θ):
    """Transform to constrained space."""
    logit_π, μ_0, log_σ_0, μ_1, log_σ_1 = θ
    π = 1. / (1. + np.exp(logit_π))
    σ_0 = np.exp(log_σ_0)
    σ_1 = np.exp(log_σ_1)
    return π, μ_0, σ_0, μ_1, σ_1

Mixture model likelihood

Suppose we observe $N$ IID realizations of $Y$, the mixture model. Denote these realizations $\boldsymbol{y} \equiv (y_1, \ldots, y_N)$.

# True parameters
π   = 0.5
μ_0 = 1.0
σ_0 = 0.7
μ_1 = 5.0
σ_1 = 1.0

N = 1_000

np.random.seed(0)
Δ   = np.random.uniform(size=N) < π
Y_0 = σ_0 * np.random.randn(N) + μ_0
Y_1 = σ_1 * np.random.randn(N) + μ_1
Y   = (1. - Δ) * Y_0 + Δ * Y_1

Given an arbitrary choice of parameters $\theta$ (not necessarily equal to the true parameters), the likelihood function gives us a way to quantify their goodness:

(we have, with a slight abuse of notation, used $f_\theta$ to denote both a joint and marginal density). For the mixture model,

where $\phi_i$ is a normal PDF with mean $\mu_i$ and scale $\sigma_i$.

Since the likelihood is a product of probabilities (numbers between zero and one), it can cause numerical issues in the presence of a large sample size $N$. As such, we work with the log-likelihood $\ell(\theta) \equiv \log L(\theta)$ instead. A function to compute $-\ell(\theta) / N$ for the mixture model is provided below. This quantity, referred to herafter as the NLL for short, is the negative of the log-likelihood normalized by the number of samples. While normalization does not materially affect the optimization, it does allow us to compare likelihoods across different population sizes. The negation is purely a matter of convention coupled to the optimization literature in which authors prefer to pose problems as minimizations, not maximizations.

def nll(θ):
    """Compute the NLL for the mixture model."""
    π, μ_0, σ_0, μ_1, σ_1 = constrain(θ)
    φ_0 = stats.norm.pdf(Y, loc=μ_0, scale=σ_0)
    φ_1 = stats.norm.pdf(Y, loc=μ_1, scale=σ_1)
    return -np.log(φ_0 * (1. - π) + φ_1 * π).mean()

Indeed, if we evalute the NLL using the same parameters that were used to generate the data, we obtain a small number. This is evidence that the parameters explain the data well:

nll(unconstrain(π, μ_0, σ_0, μ_1, σ_1))

1.8705724767802003

On the other hand, different parameters may result in a much larger NLL. This is evidence that these alternate parameters explain the data poorly:

nll(unconstrain(π, μ_0 + 1., σ_0, μ_1 + 1., σ_1))

2.540302153968016

Classical optimization methods

Since the NLL is simply a function from $\mathbb{R}^{d}$ to $\mathbb{R}$, we can use any one of a number of classical optimization methods to minimize it. For example, Newton’s method applied to the NLL results in the iterates

where $\nabla \ell$ and $\boldsymbol{H}_\ell$ are the gradient and Hessian of the log-likelihood, respectively.

In the presence of higher dimensional parameter spaces, a quasi-Newton method should be preferred as such methods do not explicitly compute and store the Hessian. An example of such a method is BFGS, used below to approximate the MLE:

θ_0 = unconstrain(0.5, Y.min(), Y.std(), Y.max(), Y.std())
result = optimize.minimize(nll, θ_0, method="BFGS")

As shown in the figure below, the convergence of BFGS is somewhat undesirable:

• the NLL is not necessarily decreased on each step and
• the convergence rate is erratic.

We cover expectation maximization, which does not suffer these ailments, in the next section.

Expectation maximization (EM)

For brevity, let $Z$ denote the latent variables. As usual, we use boldface notation to refer to multiple realizations (e.g., $\boldsymbol{y}$).

EM is a fixed point iteration on the map

Proposition. The iterates produced by EM have nondecreasing likelihood.

Proof. First, note that, by the properties of the logarithm,

Let $\theta$ be arbitrary parameter. Since $T\theta$ is a maximizer, it follows that

where the last inequality is a consequence of the nonnegativity of the Kullback-Leibler divergence. Since $L(\theta) \equiv f_\theta(\boldsymbol{y})$ by definition, the desired result follows. $\blacksquare$

Returning to the mixture model, in order to set up the EM iteration, we first have to obtain an expression for $T\theta$. For the mixture model, it will be helpful to consider only the coin flip as latent, and hence we set $Z \equiv \Delta$. We denote realizations of $\Delta$ by the lower case letter $\delta$. The joint density for a single sample $(y,\delta)$ is

Therefore,

Moreover, since

it follows that

Putting this all together,

Next, let

Note that

Therefore, differentiating $T(\theta,\theta^{\prime})$ yields

Letting $\theta^{[k]}$ denote the $k$-th iterate of EM and $\boldsymbol{\xi}^{[k]}\equiv(\xi_{1}(\theta^{[k]}),\ldots,\xi_{N}(\theta^{[k]}))^{\intercal}$, it follows that the equations for the $(k+1)$-th iterate are

def em_update(π, μ_0, σ_0, μ_1, σ_1):
    φ_0 = stats.norm.pdf(Y, loc=μ_0, scale=σ_0)
    φ_1 = stats.norm.pdf(Y, loc=μ_1, scale=σ_1)
    # The code below is optimized for readability, not speed
    ξ = φ_1 * π / (φ_0 * (1. - π) + φ_1 * π)
    π = ξ.mean()
    μ_0 = Y @ (1. - ξ) / (1. - ξ).sum()
    σ_0 = np.sqrt((1. - ξ) @ (Y - μ_0)**2 / (1. - ξ).sum())
    μ_1 = Y @ ξ / ξ.sum()
    σ_1 = np.sqrt(ξ @ (Y - μ_1)**2 / ξ.sum())
    return π, μ_0, σ_0, μ_1, σ_1

# Use the same initial guess as in Newton's method
π_k, μ_0_k, σ_0_k, μ_1_k, σ_1_k = constrain(θ_0)
for _ in range(10):
    π_k, μ_0_k, σ_0_k, μ_1_k, σ_1_k = em_update(π_k, μ_0_k, σ_0_k, μ_1_k, σ_1_k)

As shown below, EM converges in just a few iterations. Near convergence, the NLL can be non-monotone up to noise due purely to floating point error.

tabulate(
    [
        ("Exact", π, μ_0, σ_0, μ_1, σ_1),
        ("BFGS", ) + constrain(result.x),
        ("EM", π_k, μ_0_k, σ_0_k, μ_1_k, σ_1_k)
    ],
    headers=("", "π", "μ_0", "σ_0", "μ_1", "σ_1"),
    tablefmt="html",
)

Markov Chain Monte Carlo (MCMC)

So far we have treated true parameter $\theta$ as a constant. Bayesian methods treat it as a random variable with prior density $f(\theta)$. Priors express our beliefs about the parameter before seeing (i.e., unconditional on) data. After data is observed, we obtain the posterior, reflecting our updated beliefs:

That is, the posterior is proportional to the likelihood times the prior:

One can imagine any number of Bayesian estimators constructed from the posterior. Any one of the mean, median, or mode of the posterior produce reasonable estimators.

Below is an example using PyMC to compute a highest density interval (HDI), a form of credible interval, for the parameter dimensions. Reasonable priors for each dimension are chosen (e.g., $\pi \sim \operatorname{Uniform}(0, 1)$).

with pymc.Model() as model:
    π_rv = pymc.Uniform("π", lower=0., upper=1.)
    μ_0_rv = pymc.Normal("μ_0", mu=np.quantile(Y, 0.25), sigma=1.)
    σ_0_rv = pymc.HalfNormal("σ_0", sigma=1.)
    μ_1_rv = pymc.Normal("μ_1", mu=np.quantile(Y, 0.75), sigma=1.)
    σ_1_rv = pymc.HalfNormal("σ_1", sigma=1.)
    Y_0_dist = pymc.Normal.dist(mu=μ_0_rv, sigma=σ_0_rv)
    Y_1_dist = pymc.Normal.dist(mu=μ_1_rv, sigma=σ_1_rv)
    likelihood = pymc.Mixture("Y", w=[π_rv, 1. - π_rv], comp_dists=[Y_0_dist, Y_1_dist], observed=Y)
    trace = pymc.sample(2000)

pymc.plot_posterior(trace, var_names=["π", "μ_0", "σ_0", "μ_1", "σ_1"]);