{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# Prior distributions\n", "\n", "One important aspect of Bayesian inference has not yet been discussed in this tutorial: [prior distributions](https://en.wikipedia.org/wiki/Prior_probability). In Bayesian statistics, one has to provide probability (density) values for every possible parameter value *before* taking into account the data at hand. This prior distribution thus reflects all *prior* knowledge of the system that is to be investigated. In the case that no prior knowledge is available, a *non-informative* prior in the form of the so-called [Jeffreys prior](https://en.wikipedia.org/wiki/Jeffreys_prior) allows to minimize the effect of the prior on the results. The next two sub-sections discuss how one can set custom prior distributions for the parameters of the observation model and for hyper-parameters in a hyper-study or change-point study." ] }, { "cell_type": "code", "execution_count": 1, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-04-27T19:17:29.386451Z", "iopub.status.busy": "2026-04-27T19:17:29.386019Z", "iopub.status.idle": "2026-04-27T19:17:30.357735Z", "shell.execute_reply": "2026-04-27T19:17:30.357178Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Created new study.\n", "+ Successfully imported example data.\n" ] } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt # plotting\n", "plt.style.use('seaborn-v0_8-whitegrid') # plot styling\n", "\n", "import numpy as np\n", "import bayesloop as bl\n", "\n", "# prepare study for coal mining data\n", "S = bl.Study()\n", "S.load_example_data()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Parameter prior\n", "\n", "*bayesloop* employs a forward-backward algorithm that is based on [Hidden Markov models](http://www.cs.sjsu.edu/~stamp/RUA/HMM.pdf). This inference algorithm iteratively produces a parameter distribution for each time step, but it has to start these iterations from a specified probability distribution - the parameter prior. All built-in observation models already have a predefined prior, stored in the attribute `prior`. Here, the prior distribution is stored as a Python function that takes as many arguments as there are parameters in the observation model. The prior distributions can be looked up directly within `observation_models.py`. For the `Poisson` model discussed in this tutorial, the default prior distribution is defined in a method called `jeffreys` as\n", "```\n", "def jeffreys(x):\n", " return np.sqrt(1. / x)\n", "```\n", "corresponding to the non-informative Jeffreys prior, $p(\\lambda) \\propto 1/\\sqrt{\\lambda}$. *bayesloop* can also attempt to derive this type of prior for symbolic user-defined observation models; see the [custom observation model tutorial](customobservationmodels.ipynb).\n", "\n", "### Prior functions and arrays\n", "\n", "To change the predefined prior of a given observation model, one can add the keyword argument `prior` when defining an observation model. There are different ways of defining a parameter prior in *bayesloop*: If `prior=None` is set, *bayesloop* will assign equal probability to all parameter values, resulting in a uniform prior distribution within the specified parameter boundaries. One can also directly supply a NumPy array with prior probability (density) values. The shape of the array must match the shape of the parameter grid! Another way to define a custom prior is to provide a function that takes exactly as many arguments as there are parameters in the defined observation model. *bayesloop* will then evaluate the function for all parameter values and assign the corresponding probability values.\n", "\n", "
\n", "**Note:** In all of the cases described above, *bayesloop* will re-normalize the provided prior values, so they do not need to be passed in a normalized form. Below, we describe the possibility of using probability distributions from the SymPy stats module as prior distributions, which are not re-normalized by *bayesloop*.\n", "
\n", "\n", "Next, we illustrate the difference between the Jeffreys prior and a flat, uniform prior with a very simple inference example: We fit the coal mining example data set using the `Poisson` observation model and further assume the rate parameter to be static:" ] }, { "cell_type": "code", "execution_count": 2, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-04-27T19:17:30.359138Z", "iopub.status.busy": "2026-04-27T19:17:30.358954Z", "iopub.status.idle": "2026-04-27T19:17:30.437176Z", "shell.execute_reply": "2026-04-27T19:17:30.436762Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Transition model: Static/constant parameter values. Hyper-Parameter(s): []\n", "Fit with built-in Jeffreys prior:\n", "+ Observation model: Poisson. Parameter(s): ['accident_rate']\n", "+ Started new fit:\n", " + Formatted data.\n", " + Set prior (function): jeffreys. Values have been re-normalized.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "2f21430dd37e4877ab2d02ffec6c2135", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/110 [00:00. Values have been re-normalized.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "19ff71884d594379acc0e02ebcc8d40f", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/110 [00:00\n", "**Note:** The support interval of a prior distribution defined via SymPy can deviate from the parameter interval specified in *bayesloop*. In the example above, we specified the parameter interval ]0, 6[, while the exponential prior has the support ]0, $\\infty$[. SymPy priors are not re-normalized with respect to the specified parameter interval. Be aware that the resulting model evidence value will only be correct if no parameter values outside of the parameter boundaries gain significant probability values. In most cases, one can simply check whether the parameter distribution has sufficiently *fallen off* at the parameter boundaries.\n", "\n", "\n", "## Hyper-parameter priors\n", "\n", "As shown before, [hyper-studies](hyperstudy.ipynb) and [change-point studies](changepointstudy.ipynb) can be used to determine the full distribution of hyper-parameters (the parameters of the transition model). As for the time-varying parameters of the observation model, one might have prior knowledge about the values of certain hyper-parameters that can be included into the study to refine the resulting distribution of these hyper-parameters. Hyper-parameter priors can be defined just as regular priors, either by an arbitrary function or by a list of `sympy.stats` random variables.\n", "\n", "In a first example, we return to the simple change-point model of the coal-mining data set and perform two fits of the change-point: first, we specify no hyper-prior for the time step of our change-point, assuming equal probability for each year in our data set. Second, we define a Normal distribution around the year 1920 with a (rather unrealistic) standard deviation of 5 years as the hyper-prior using a SymPy random variable. For both fits, we plot the change-point distribution to show the differences induced by the different priors:" ] }, { "cell_type": "code", "execution_count": 5, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-04-27T19:17:30.460216Z", "iopub.status.busy": "2026-04-27T19:17:30.460120Z", "iopub.status.idle": "2026-04-27T19:17:31.268419Z", "shell.execute_reply": "2026-04-27T19:17:31.267957Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "Fit with flat hyper-prior:\n", "+ Created new study.\n", " --> Hyper-study\n", " --> Change-point analysis\n", "+ Successfully imported example data.\n", "+ Observation model: Poisson. Parameter(s): ['accident_rate']\n", "+ Transition model: Change-point. Hyper-Parameter(s): ['tChange']\n", "+ Detected 1 change-point(s) in transition model: ['tChange']\n", "+ Set hyper-prior(s): ['uniform']\n", "+ Started new fit.\n", " + 109 analyses to run.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "02dd3926b5204860955b0c844e608306", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/109 [00:00" ] }, "metadata": {}, "output_type": "display_data" }, { "name": "stdout", "output_type": "stream", "text": [ "-----\n", "\n", "Fit with custom normal prior:\n", "+ Transition model: Change-point. Hyper-Parameter(s): ['tChange']\n", "+ Detected 1 change-point(s) in transition model: ['tChange']\n", "+ Set hyper-prior(s): ['sqrt(2)*exp(-(x - 1920)**2/50)/(10*sqrt(pi))']\n", "+ Started new fit.\n", " + 109 analyses to run.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "b8ec5217fb6c47509e414c9398f8aa35", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/109 [00:00" ] }, "metadata": {}, "output_type": "display_data" } ], "source": [ "print('Fit with flat hyper-prior:')\n", "S = bl.ChangepointStudy()\n", "S.load_example_data()\n", "\n", "L = bl.om.Poisson('accident_rate', bl.oint(0, 6, 1000))\n", "T = bl.tm.ChangePoint('t_change', 'all')\n", "\n", "S.set(L, T)\n", "S.fit()\n", "\n", "plt.figure(figsize=(8,4))\n", "S.plot('t_change', facecolor='g', alpha=0.7)\n", "plt.xlim([1870, 1930])\n", "plt.show()\n", "print('-----\\n')\n", " \n", "print('Fit with custom normal prior:')\n", "T = bl.tm.ChangePoint('t_change', 'all', prior=sympy.stats.Normal('norm', 1920, 5))\n", "S.set(T)\n", "S.fit()\n", "\n", "plt.figure(figsize=(8,4))\n", "S.plot('t_change', facecolor='g', alpha=0.7)\n", "plt.xlim([1870, 1930]);" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Since we used a quite narrow prior (containing a lot of information) in the second case, the resulting distribution is strongly shifted towards the prior. The following example revisits the two-break-point model from [here](changepointstudy.ipynb#Analyzing-structural-breaks-in-time-series-models) and a linear decrease with a varying slope as a hyper-parameter. Here, we define a Gaussian prior for the slope hyper-parameter, which is centered around the value -0.2 with a standard deviation of 0.4, via a lambda-function. For simplification, we set the break-points to fixed years." ] }, { "cell_type": "code", "execution_count": 6, "metadata": { "collapsed": false, "execution": { "iopub.execute_input": "2026-04-27T19:17:31.269755Z", "iopub.status.busy": "2026-04-27T19:17:31.269642Z", "iopub.status.idle": "2026-04-27T19:17:31.481833Z", "shell.execute_reply": "2026-04-27T19:17:31.481359Z" } }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "+ Created new study.\n", " --> Hyper-study\n", "+ Successfully imported example data.\n", "+ Observation model: Poisson. Parameter(s): ['accident_rate']\n", "+ Transition model: Serial transition model. Hyper-Parameter(s): ['slope', 't_1', 't_2']\n", "+ Set hyper-prior(s): [' (re-normalized)', 'uniform', 'uniform']\n", "+ Started new fit.\n", " + 30 analyses to run.\n" ] }, { "data": { "application/vnd.jupyter.widget-view+json": { "model_id": "710fcd390f9b4f55a92ff78e29a9081e", "version_major": 2, "version_minor": 0 }, "text/plain": [ " 0%| | 0/30 [00:00