Commit 1cc5fc49 authored by Gregory Ashton's avatar Gregory Ashton

Clean up of example and add explanation

parent 7bf26c23
Pipeline #22768 passed with stages
in 11 minutes and 11 seconds
As part of the :code:`tupak.result.Result` object, we provide a method to
calculate the Occam factor (c.f., Chapter 28, `Mackay "Information Theory,
Inference, and Learning Algorithms"
<>`). This is an approximate
estimate based on the posterior samples, and assumes the posteriors are well
approximate by a Gaussian.
The Occam factor penalizes models with larger numbers of parameters (or
equivalently a larger "prior volume"). This example won't try to go through
explaining the meaning of this, or how it is calculated (those details are
sufficiently well done in Mackay's book linked above). Insetad, it demonstrates
how to calculate the Occam factor in :code:`tupak` and shows an example of it
working in practise.
If you have a :code:`result` object, the Occam factor can be calculated simply
from :code:`result.occam_factor(priors)` where :code:`priors` is the dictionary
of priors used during the model fitting. These priors should be uniform
priors only. Other priors may cause unexpected behaviour.
In the example, we generate a data set which contains Gaussian noise added to a
quadratic function. We then fit polynomials of differing degree. The final plot
shows that the largest evidence favours the quadratic polynomial (as expected)
and as the degree of polynomial increases, the evidence falls of in line with
the increasing (negative) Occam factor.
Note - the code uses a course 100-point estimation for speed, results can be
improved by increasing this to say 500 or 1000.
from __future__ import division
import tupak
......@@ -70,24 +98,32 @@ def fit(n):
result = tupak.run_sampler(
likelihood=likelihood, priors=priors, npoints=100, outdir=outdir,
return result.log_evidence, result.log_evidence_err, np.log(result.occam_factor(priors))
return (result.log_evidence, result.log_evidence_err,
fig, ax = plt.subplots()
fig, ax1 = plt.subplots()
ax2 = ax1.twinx()
log_evidences = []
log_evidences_err = []
log_occam_factors = []
ns = range(1, 10)
ns = range(1, 11)
for l in ns:
e, e_err, o = fit(l)
ax.errorbar(ns, log_evidences-np.max(log_evidences), yerr=log_evidences_err,
fmt='-o', color='C0')
ax.plot(ns, log_occam_factors-np.max(log_occam_factors),
'-o', color='C1', alpha=0.5)
ax1.errorbar(ns, log_evidences, yerr=log_evidences_err,
fmt='-o', color='C0')
ax1.set_ylabel("Unnormalized log evidence", color='C0')
ax1.tick_params('y', colors='C0')
ax2.plot(ns, log_occam_factors,
'-o', color='C1', alpha=0.5)
ax2.tick_params('y', colors='C1')
ax2.set_ylabel('Occam factor', color='C1')
ax1.set_xlabel('Degree of polynomial')
fig.savefig('{}/{}_test'.format(outdir, label))
......@@ -210,8 +210,9 @@ class Result(dict):
def occam_factor(self, priors):
""" The Occam factor,
See Mackay "Information Theor, Inference and Learning Algorithms,
Cambridge University Press (2003) Eq. (28.10)
See Chapter 28, `Mackay "Information Theory, Inference, and Learning
Algorithms" <>`_ Cambridge
University Press (2003).
return self.posterior_volume / self.prior_volume(priors)
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment