I keep meaning to write a post on this...

At least within Bayesian probability, there is a single unique unambiguously-correct answer to "how should we penalize for complexity?": the Bayes factor. The Bayes factor is Hard to compute in general, which is why there's a whole slew of of other numbers which approximate it in various ways.

Here's how Bayes factors work. Want to know whether model 1 or model 2 is more consistent with the data? Then compute $P\left[mode{l}_1|data\right]$ and $P\left[mode{l}_2|data\right]$. Using Bayes' rule:

$P\left[mode{l}_i|data\right]=\frac{1}{Z}P\left[data|mode{l}_i\right]P\left[mode{l}_i\right]$

where Z is the normalizer. If we're just comparing two models, then we can get rid of that annoying Z by computing odds for the two models:

$\frac{P\left[mode{l}_1|data\right]}{P\left[mode{l}_2|data\right]}=\frac{P\left[data|mode{l}_1\right]}{P\left[data|mode{l}_2\right]}\frac{P\left[mode{l}_1\right]}{P\left[mode{l}_2\right]}$

In English: posterior relative odds of the two models is equal to prior odds times the ratio of likelihoods. That likelihood ratio $\frac{P\left[data|mode{l}_1\right]}{P\left[data|mode{l}_2\right]}$ is the Bayes' factor: it directly describes the update in the relative odds of the two models, due to the data.

Example

20 coin flips yield 16 heads and 4 tails. Is the coin biased?

Here we have two models:

• Model 1: coin unbiased
• Model 2: coin has some unknown probability $\theta$ of coming up heads (we'll use a uniform prior on $\theta$ for simplicity; this is not a good idea in general)

The second model has one free parameter (the bias) which we can use to fit the data, but it’s more complex and prone to over-fitting. When we integrate over that free parameter, it will fit the data poorly over most of the parameter space - thus the "penalty" associated with free parameters in general.

In this example, the integral is exactly tractable (it's a dirichlet-multinomial model), and we get:

• $P\left[data|mode{l}_1\right]=(\frac{20}{4})\left(\frac{1}{2}{)}^{16}\right(\frac{1}{2}{)}^4=0.0046$
• $P\left[data|mode{l}_2\right]={\int }_{\theta }(\frac{20}{4})\left(\theta {)}^{16}\right(1-\theta {)}^{4}d\theta =0.048$

So the Bayes factor is (.048)/(.0046) ~ 10, in favor of a biased coin. In practice, I'd say unbiased coins are at least 10x more likely than biased coins in day-to-day life a priori, so we might still think the coin is unbiased. But if we were genuinely unsure to start with, then this would be pretty decent evidence in favor.

Approximation

Why is this hard to compute in general? Well, if the model has an m-dimensional vector of free parameters $\theta$, then $P\left[data|model\right]={\int }_{\theta }P[data|model,\theta ]dP\left[\theta |model\right]$, which is an integral in m dimensions. That's not just NP-complete, it's #P-complete.

There's various ways to approximate that integral. The two most common are:

• MCMC (see e.g. here for application to Bayes factors): main limitation is that MCMC gets finicky in high dimensions.
• Laplace approximation around the maximum-likelihood point. Aside from the usual limitations of using a single maximum-likelihood point, main limitation is that we need the determinant of the Hessian, which takes somewhere between O(k) and O($k^3$) to compute for most k-free-parameter models (depending on how clever we are with the linear algebra).

BIC

The BIC can be derived directly from the Laplace approximation. Laplace says:

${\int }_{\theta }P\left[data|\theta \right]p\left[\theta \right]d\theta \approx P\left[data|{\theta }^{\ast }\right]p\left[{\theta }^{\ast }\right]\left(2\pi {)}^{\frac{k}{2}}det\right(-\frac{{\partial }^{2}lnP\left[data|\theta \right]}{\partial {\theta }^2}{)}^{-\frac{1}{2}}$

where ${\theta }^{\ast }$ is the maximum likelihood point and k is the dimension of $\theta$.

Notice that $ln\left(P\right[data|{\theta }^{\ast }\left]p\right[{\theta }^{\ast }\left]\right)$ is the usual "score" of a model; the rest is what we can think of as a complexity penalty, and that's what the BIC will approximate. The main assumption is that $\frac{{\partial }^{2}lnP\left[data|\theta \right]}{\partial {\theta }^2}$ scales with the number of data points N; this is trivial if the data points are iiid (given $\theta$), since in that case $lnP\left[data|\theta \right]$ is a sum over data points.

Now for the trick: that's the only part of the Laplace approximation's "complexity penalty" which scales with the number of data points N. It scales in proportion to $N^{-\frac{k}{2}}$, since each of k dimensions of the determinant contributes a factor N. Take the log of the whole thing, and all the rest will be constants which are dominated by the $\frac{k}{2}lnN$ term as N goes to infinity. And that's the BIC: log of the MAP probability, minus $\frac{k}{2}lnN$.

Finally, the important part: the main benefit of knowing all this is being able to predict when the BIC will and will not work well. If it's an approximation of Laplace, which is itself an approximation of the Bayes factor, then BIC will work well exactly when those approximations are accurate. In particular, BIC fails when the terms we ignored (i.e. $ln\left(\right(2\pi {)}^{\frac{k}{2}}det(-\frac{1}{N}\frac{{\partial }^{2}lnP\left[data|\theta \right]}{\partial {\theta }^2}{)}^{-\frac{1}{2}})$) are as large as or even larger than $ln\left(P\right[data|{\theta }^{\ast }\left]p\right[{\theta }^{\ast }\left]\right)-\frac{k}{2}lnN$, i.e. the BIC itself.

A different view is to look at the search process for the models, rather than the model itself. If model A is found from a process that evaluates 10 models, and model B is found from a process that evaluates 10,000, and they otherwise have similar results, then A is much more likely to generalize to new data points than B.

The formalization of this concept is called VC dimension and is a big part of Machine Learning Theory (although arguably it hasn't been very helpful in practice): https://en.wikipedia.org/wiki/Vapnik%E2%80%93Chervonenkis_dimension

The exact Bayesian solution penalizes complex models as a side effect. Each model should have a prior over its parameters. The more complex model can fit the data better, so P(data | best-fit parameters, model) is higher. But the model gets penalized because P(best-fit parameters | model) is lower on the prior. Why? The prior is thinly spread over a higher dimensional parameter space so it is lower for any particular set of parameters. This is called "Bayesian Occam's razor".

If one has 2 possible models to fit a data set, by how much should one penalize the model which has an additional free parameter?

Penalization might not be necessary if your learning procedure is stochastic and favors simple explanations. I encourage you to take a look on the nice poster/paper «Deep learning generalizes because the parameter-function map is biased towards simple functions» (PAC-Bayesian learning theory + empirical intuitions).