report.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from IPython.display import display, HTML\n",
    "import os"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Bayesian Change-point Detection on Coal Mining Disasters"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Introduction\n",
    "\n",
    "The study of time series, i.e. sequences of data points ordered in time describing the behavior of systems, dates back to the work of G. U. Yule in the 1920s [1]. The properties of such sequences can change over time due to external events and/or internal systematic changes in the underlying system dynamics or inherent distribution [2]. Therefore, a natural way of improving our understanding of the underlying systems is to partition these sequences into smaller sub-sequences, each capturing a distinct set of properties. The boundaries between successive partitions are known as change-points, and the problem of identifying their occurrence and pinpointing their location in time is known as change-point detection [3]. \n",
    "\n",
    "There is always a fine balance between misidentifying a regular fluctuation as a change in properties, and overlooking a significant change when it occurs. This balance becomes even harder to maintain in scenarios where the factors influencing the data are unknown, which is exactly where Bayesian methods shine [4]. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots/data_plot.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(f\"<div style='text-align: center;'><img src='plots/data_plot.png'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Problem description\n",
    "\n",
    "In this project we will perform Bayesian change-point detection on a historical data set detailing coal mining disasters in the United Kingdom from 1851 to 1962. These records, originally given by Maguire, Pearson and Wynn in 1952, and subsequently corrected and extended by Jarrett in 1979 [5], can be traced back to the Colliery Year Book and Coal Trades Directory, and have become a well-established and widely used dataset in the field of change-point analysis.\n",
    "\n",
    "The cumulative number of the said coal mining disasters is depicted in Figure 1. A key observation from this figure is that the rate of disasters did not remain constant over the 112 years. The most prominent change was observed in the late 1880s, marking a period where a previously stable rate of disasters declined considerably. This change has been attributed to certain improvements in safety regulations during that period [6], which is an interesting topic in itself. However, this project will strictly focus on determining the number of significant change-points that occurred during the 112-year period and estimating when these are most likely to have occurred, deliberately excluding investigation of the social or economic aspects behind them.\n",
    "\n",
    "We formulate our statistical model as follows.\n",
    "\n",
    "Let $ y_{1:N} := (y_1, y_2, \\ldots, y_N) $ be the series of the count of annual mining disasters, where $N = 112$ and $\\forall y_i \\in \\mathbb{R}, i \\in \\mathbb{N}$. Suppose that $ \\{S^{(1)}, S^{(2)}, \\ldots, S^{(l)}\\} $ are disjoint segments of $ y_{1:N} $ , where the number of change-points $ l \\in \\mathbb{N} $ is unknown. Define a change-point $ T_j $ to be the index of the first element of subset $ S^{(j)} $, $\\forall j \\in \\{1, \\ldots, l\\}$. The goal is to determine each set $ S^{(j)} $, and therefore the values of each $ T_j $, subject to the segmentation $ S^{(1)}, S^{(2)}, \\ldots, S^{(l)} $ being optimal according to the marginal likelihood of the model.\n",
    "\n",
    "We will consider the count of annual mining disasters $y_i$ as an observation originating from a Poisson process with an unknown rate $\\lambda_j$:\n",
    "\n",
    "$$\n",
    "y_{1:N} \\sim \n",
    "\\begin{cases} \n",
    "  \\text{Pois}(\\lambda_1), & i = 1, \\ldots, (T_1 - 1) \\\\\n",
    "  \\text{Pois}(\\lambda_2), & i = T_1, \\ldots, (T_2 - 1) \\\\\n",
    "  \\vdots & \\vdots \\\\\n",
    "  \\text{Pois}(\\lambda_{l+1}), & i = T_l, \\ldots, N, \\text{ where } N=112\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "We will use a Gamma(2, 1) prior for the unknown rates $ \\lambda_j $, and an ordered discrete uniform prior for the unknown change-points $ T_j $:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "\\lambda_j &\\sim \\text{Gamma}(2, 1) \\\\\n",
    "T_j &\\sim \\text{Uniform}(1, N), \\text{ where } N=112\n",
    "\\end{align*}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Analytical solution\n",
    "\n",
    "##### Step 1: Likelihood Function\n",
    "\n",
    "The likelihood of observing data $ y_{1:N} $ given parameters $ \\lambda $ and $ T $ is found by the product of Poisson probabilities across each segment:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    L(y \\mid \\lambda, T) &= L(y_{1:T_1-1} \\mid \\lambda_1) \\cdots L(y_{T_l:N} \\mid \\lambda_{l+1}) \\\\\n",
    "    &= \\prod_{i=1}^{T_1 - 1} \\text{Pois}(y_i \\mid \\lambda_1) \\cdots \\prod_{i=T_l}^{N} \\text{Pois}(y_i \\mid \\lambda_{l+1}) \\\\\n",
    "    &= \\prod_{i=1}^{T_1-1} \\frac{{e^{-\\lambda_1} \\lambda_1^{y_i}}}{{y_i!}} \\cdots \\prod_{i=T_l}^{N} \\frac{{e^{-\\lambda_{l+1}} \\lambda_{l+1}^{y_i}}}{{y_i!}}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "##### Step 2: Marginal Posteriors of $ T $\n",
    "\n",
    "The marginal posterior of $ T $ is found by integrating out the $ \\lambda $ parameters:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    p(T \\mid y) \n",
    "    &\\propto \\int_{\\lambda_1=0}^{\\infty} \\cdots \\int_{\\lambda_{l+1}=0}^{\\infty} L(y \\mid \\lambda, T) \\, p(\\lambda_1) \\cdots p(\\lambda_{l+1}) \\,  d\\lambda_1  \\ldots d\\lambda_{l+1} \\\\\n",
    "    &\\propto \\int_{\\lambda_1=0}^{\\infty} \\cdots \\int_{\\lambda_{l+1}=0}^{\\infty} \\left( \\prod_{i=1}^{T_1 - 1} \\frac{e^{-\\lambda_1} \\lambda_1^{y_i}}{y_i!} \\right) \\cdots \\left( \\prod_{i=T_l}^{N} \\frac{e^{-\\lambda_{l+1}} \\lambda_{l+1}^{y_i}}{y_i!} \\right) \\frac{\\lambda_1^{\\alpha-1} e^{-\\beta \\lambda_1}}{\\Gamma(\\alpha)} \\cdots \\frac{\\lambda_{l+1}^{\\alpha-1} e^{-\\beta \\lambda_{l+1}}}{\\Gamma(\\alpha)} d\\lambda_1  \\ldots d\\lambda_{l+1} \\\\\n",
    "    &\\propto \\frac{1}{\\Gamma(\\alpha)^N} \\left( \\prod_{i=1}^{T_1 - 1} \\frac{1}{y_i!} \\right) \\cdots \\left( \\prod_{i=T_l}^{N} \\frac{1}{y_i!} \\right) \\times \\\\ \n",
    "    &\\qquad \\int_{\\lambda_1=0}^{\\infty} e^{-\\lambda_1(T_1 - 1 + \\beta)} \\lambda_1^{\\left(\\sum_{i=1}^{T_1 - 1} y_i\\right) + \\alpha - 1} d\\lambda_1 \\cdots \\int_{\\lambda_{l+1}=0}^{\\infty} e^{-\\lambda_{l+1}(N - T_l + \\beta)} \\lambda_{l+1}^{\\left(\\sum_{i=T_l}^{N} y_i\\right) + \\alpha - 1} d\\lambda_{l+1} \\\\\n",
    "    &\\propto \\left( \\prod_{i=1}^{T_1 - 1} \\frac{1}{y_i!} \\right) \\cdots \\left( \\prod_{i=T_l}^{N} \\frac{1}{y_i!} \\right) \\frac{\\Gamma\\left(\\sum_{i=1}^{T_1 - 1} y_i + \\alpha\\right)}{(T_1 - 1 + \\beta)^{\\left(\\sum_{i=1}^{T_1 - 1} y_i\\right) + \\alpha}} \\cdots \\frac{\\Gamma\\left(\\sum_{i=T_l}^{N} y_i + \\alpha\\right)}{(N - T_l + \\beta)^{\\left(\\sum_{i=T_l}^{N} y_i\\right) + \\alpha}} \\\\\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where the parameters of the prior distribution of each $\\lambda$ will be considered to be $ a = 2 $ and $ \\beta = 1 $.\n",
    "\n",
    "However, due to the non-conjugacy introduced by $ T $, this integral does not have a closed form and must be evaluated numerically for each value of $ T $.\n",
    "\n",
    "##### Step 3: Conditional Posteriors of $ \\lambda $\n",
    "\n",
    "For each fixed (sampled) value $ T_j $, we can calculate the conditional posterior of $ \\lambda_j $ by considering only the data points that fall into the corresponding segment $ S^{(j)} $:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    p(\\lambda_j \\mid y_{T_{j-1}:T_j-1}) \n",
    "    &\\propto L(y_{T_{j-1}:T_j-1} \\mid \\lambda_j) \\, p(\\lambda_j) \\quad \\text{for } j = 1, \\ldots N \\\\\n",
    "    &\\propto \\left( \\prod_{i=T_{j-1}}^{T_j-1} e^{-\\lambda_j} \\lambda_j^{y_i} \\right) \\beta^\\alpha \\lambda_j^{\\alpha-1} e^{-\\beta \\lambda_j} \\\\\n",
    "    &\\propto e^{-\\lambda_j(T_j-1-T_{j-1}+\\beta)} \\lambda_j^{\\sum_{i=T_{j-1}}^{T_j-1} y_i+\\alpha-1}\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "Recognizing that this is the kernel of a Gamma distribution, we can sample parameter $ \\lambda_j $ as:\n",
    "\n",
    "$$\n",
    "\\lambda_j \\mid y_{T_{j-1}:T_j-1} \\sim Gamma\\left(\\sum_{i=T_{j-1}}^{T_j-1} y_i + \\alpha, \\, T_j-1-T_{j-1}+\\beta \\right)\n",
    "$$\n",
    "\n",
    "This is expected because the Gamma distribution is a conjugate prior for the Poisson distribution, and therefore the posterior distribution for $ \\lambda_j $ should also be a Gamma distribution, with the initial Gamma parameters updated taking into account the relevant data: \n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "    \\alpha_{new} &= \\alpha + \\sum_{i \\in S^{(j)}} y_i \\\\\n",
    "    \\beta_{new} &= \\beta + L_j\n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "where $ Q_j = \\sum_{i \\in S^{(j)}} y_i = \\sum_{i=T_{j-1}}^{T_j-1} y_i $ represents the sum of the data-points, i.e. the disaster counts, within segment $ S^{(j)} $, and $ L_j = T_j-1-T_{j-1} $ represents the length of $ S^{(j)} $.\n",
    "\n",
    "##### Step 4: Posterior Predictive Distribution\n",
    "\n",
    "Using the sampled values of $ T $ and $ \\lambda $, we can then generate the posterior predictive distribution for new data points.\n",
    "\n",
    "This entire process is computationally intensive and needs to be executed using a numerical method, in our case MCMC.\n",
    "\n",
    "##### Step 5: Marginal Likelihood\n",
    "\n",
    "The marginal likelihood of a model with $l$ change-points can be evaluated as:\n",
    "\n",
    "$$\n",
    "\\begin{align*}\n",
    "p(M \\mid y) \n",
    "&= \\sum_{T_1=2}^{N+1-l} \\sum_{T_2=T_1+1}^{N+2-l} \\cdots \\sum_{T_l=T_{l-1}+1}^{N} \\int_{\\lambda_1=0}^{\\infty} \\cdots \\int_{\\lambda_{l+1}=0}^{\\infty} p(y \\mid T_1,\\ldots,T_l,\\lambda_1,\\ldots,\\lambda_{l+1}) p(T_1,\\ldots,T_l,\\lambda_1,\\ldots,\\lambda_{l+1}) \\,d\\lambda_1 \\ldots d\\lambda_{l+1} \\\\\n",
    "&= \\sum_{T_1=2}^{N+1-l} \\sum_{T_2=T_1+1}^{N+2-l} \\cdots \\sum_{T_l=T_{l-1}+1}^{N} \\binom{N-1}{l}^{-1} \\frac{\\prod_{j=1}^{l+1} (L_j + 1)^{-Q_j-2} \\Gamma (2 + Q_j)}{\\prod_{i=1}^{N} y_i!} \n",
    "\\end{align*}\n",
    "$$\n",
    "\n",
    "\n",
    "These nested sums are computationally intensive for a large number of change-points $l$, which again begs the need to employ a numerical method for their calculation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### MCMC implementation\n",
    "\n",
    "The intractable posterior distributions (described in the previous section) were approximated using PyMC's [7] implementation of Sequential Monte Carlo sampling (SMC).  \n",
    "\n",
    "SMC methods are a general class of Monte Carlo methods tailored to solve statistical (predominantly Bayesian) inference problems recursively. Specifically, they approximate the posterior distribution of a probabilistic model by sampling iteratively from a sequence of target probability densities [8], [9]. To avoid the problem of sampling from difficult target probability density functions, standard MCMC samplers can be converted into SMC samplers, as was presented for the first time in [9] and [10]. Specifically, PyMC bases its SMC implementation on CATMIP [11], which combines tempering (i.e. iterating through a sequence of target distributions) and resampling, with the standard MCMC simulation of the Metropolis algorithm [12].\n",
    "\n",
    "**Algorithm Summary**\n",
    "\n",
    "SMC sampling works by moving through successive stages. At each stage, samples are generated from a tempered posterior: \n",
    "$$ p(\\theta \\mid y)_\\beta = p(y \\mid \\theta)^\\beta p(\\theta) $$\n",
    "This is a \"tempered\" form of the posterior distribution $ p(\\theta | y) $, where $ \\beta $ (the tempering parameter) moves from 0 (prior distribution) to 1 (true posterior distribution), according to a certain adaptive process.\n",
    "\n",
    "- **Step 1:** Initialize $ \\beta $ at zero and start at stage zero.\n",
    "- **Step 2:** Generate $N$ samples $S_\\beta$ from the prior (when $\\beta = 0$ the tempered posterior is the prior).\n",
    "- **Step 3:** Increase $\\beta$ in order to make the effective sample size equal to some predefined value (typically half the total number of samples). This is to allow for a gradual increase in the influence of the likelihood on the posterior distribution without losing the diversity of the sample set, which is critical for an accurate approximation of the posterior. \n",
    "The process of updating $ \\beta $ typically involves the following steps:\n",
    "   1. **Increment $ \\beta $**: Increase $ \\beta $ by a small amount to $ \\beta_{\\text{new}} $.\n",
    "   2. **Calculate Weights for $ \\beta_{\\text{new}} $**:\n",
    "      Compute the raw weights $ W_i $ for each sample for the current $ \\beta_{\\text{new}} $, based on the following likelihood ratio:\n",
    "      $$ \n",
    "         W_i = \\frac{p(y | \\theta_i)^{\\beta_{\\text{new}}} p(\\theta_i)}{p(y | \\theta_i)^{\\beta_{\\text{old}}} p(\\theta_i)} = p(y | \\theta_i) ^{\\beta_{\\text{new}} - \\beta_{\\text{old}}} \n",
    "      $$\n",
    "   3. **Normalize the Weights**:\n",
    "      Normalize the raw weights so that they sum up to 1:\n",
    "      $$ \n",
    "         w_i = \\frac{w_i(\\beta)}{\\sum_{j=1}^{N} w_j(\\beta)} \n",
    "      $$\n",
    "   4. **Calculate the ESS**:\n",
    "      With these normalized weights, calculate the ESS as:\n",
    "      $$ \n",
    "         \\text{ESS} = \\frac{\\left( \\sum_{i=1}^{N} w_i \\right)^2}{\\sum_{i=1}^{N} w_i^2} \n",
    "      $$\n",
    "   5. **Check ESS and repeat if possible**:\n",
    "       If the ESS at $ \\beta_{\\text{new}} $ is above the threshold, use $ \\beta_{\\text{new}} $ as the new $ \\beta $ value for the next iteration of the SMC algorithm. Otherwise, go back to step 1 to calculate a new $ \\beta_{\\text{new}} $.\n",
    "- **Step 4:** Obtain a new set of samples $ S_{w} $ by re-sampling $ S_{\\beta} $ according to $ W $ ($ W $ is essentially the last set of weights that was computed in the previous step).\n",
    "- **Step 5:** Use $ S_{w} $ to compute the mean and covariance for the proposal distribution, which is typically assumed to be multivariate normal. This step is crucial because it adjusts the proposal distribution to better match the shape of the target distribution. By using the covariance of the samples, the proposal distribution can make more informed proposals, leading to higher acceptance rates in the next steps. \n",
    "- **Step 6:** Run $N$ independent MCMC (in our case Metropolis-Hastings) chains, starting each one from a different sample in $S_w$. \n",
    "   - In the Metropolis-Hastings algorithm, each chain proposes a new sample $ \\theta^* $ based on a proposal distribution $ q(\\theta^*|\\theta_i) $, often a multivariate normal distribution centred on the current sample with covariance $ c \\times \\Sigma $:\n",
    "   $$ \\theta^* \\sim \\mathcal{N}(\\theta_i, c \\times \\Sigma) $$\n",
    "   - The proposed sample is accepted with the probability:\n",
    "   $$ \\alpha(\\theta_i, \\theta^*) = \\text{min} \\left(1, \\frac{p(\\theta^*|y)q(\\theta_i|\\theta^*)}{p(\\theta_i|y)q(\\theta^*|\\theta_i)} \\right) $$   \n",
    "   - If a uniform random number $ u \\sim U(0,1) $ is less than $ \\alpha(\\theta_i, \\theta^*) $, the new sample $ \\theta^* $ is accepted, and the new state of the chain becomes $ \\theta^* $; otherwise, the chain remains at $ \\theta_i $.\n",
    "   - This process is iterated $ n_{steps} $ times to form a chain of samples, and the chains may be run in parallel.\n",
    "   - For stages other than 0, the acceptance rate $ \\alpha $ from the previous stage is used to adjust the scaling of the proposal distribution. If the acceptance rate is too low, it suggests the chain is moving too far from the current sample, so $ c $ is reduced to make smaller moves. Conversely, if the acceptance rate is too high, it could mean the chain is not exploring an adequate part of space, so $ c $ is increased. The new scaling factor is calculated as: $ c_{\\text{new}} = c_{\\text{old}} \\times \\text{min} \\left( \\frac{\\alpha}{\\alpha_{\\text{target}}}, f \\right) $ , where $ \\alpha_{\\text{target}} $ is the desired acceptance rate (typically around 0.23 for a single parameter, or 0.44 for high-dimensional problems) and $ f $ is a limiting factor to prevent too drastic changes in scaling.\n",
    "- **Step 7:** The $N$ chains are run until the autocorrelation with the samples from the previous stage stops decreasing given a certain threshold.\n",
    "- **Step 8:** Repeat from step 3 until $ |\\beta| \\geq 1 $.\n",
    "- **Step 9:** Collect the final set of $N$ samples from the posterior.\n",
    "\n",
    "In our case, we ran 8 independent chains for each model. To be able to achieve convergence for each model, we fine-tuned the autocorrelation threshold, i.e. the hyper-parameter that controls how long each Metropolis-Hastings chain will run. The values we ended up using for each model case are: $10^{-2}$ for 0 and 1 change-points, $ 5 \\times 10^{-3} $ for 2 change-points, $ 10^{-4} $ for 3 change-points, and $ 10^{-6} $ for 4 change-points. This fine-tuning was crucial since, as the complexity of the model increases, the algorithm needs to search the parameter space for longer to be able to adequately approximate the optimal solution.\n",
    "\n",
    "**High-dimensional setting**\n",
    "\n",
    "When dealing with many parameters each with a different prior in Sequential Monte Carlo (SMC) methods, the complexity increases but the core algorithm remains the same:\n",
    "- **Joint Prior Distribution**: Each parameter $ \\theta_i $ has a distinct prior $ p(\\theta_i) $. The joint prior distribution for all parameters is the product of these individual priors, assuming parameter independence:\n",
    "   $$ p(\\theta_1, \\theta_2, ..., \\theta_n) = \\prod_{i=1}^{n} p(\\theta_i) $$\n",
    "- **Tempering Process**: The tempering parameter $ \\beta $ incrementally shifts the focus from the joint prior to the joint posterior. This is done by tempering the likelihood:\n",
    "   $$ p(y|\\theta_1, \\theta_2, ..., \\theta_n)^\\beta $$\n",
    "   where $ y $ is the data, and $ \\theta_1, \\theta_2, ..., \\theta_n $ are the parameters.\n",
    "- **Sampling and Proposal Distribution**: Sampling is performed in the multi-dimensional space of all parameters. The proposal distribution, typically a multivariate normal, must consider the covariance between parameters. This distribution might have a complex structure, especially if parameters are correlated.\n",
    "- **Metropolis-Hastings Steps**: Each Metropolis chain operates in this multi-dimensional space. Proposals are generated for all parameters simultaneously, and their acceptance is based on the joint likelihood and the priors of all parameters.\n",
    "- **Importance Weights and Resampling**: The weights are calculated based on the full model's likelihood and priors. These weights are then used in the resampling step to focus on more probable regions of the parameter space.\n",
    "- **Computational Considerations**: The process becomes computationally more intensive with more parameters, especially if the parameter space is high-dimensional or if parameters are highly correlated.\n",
    "\n",
    "**Marginal likelihood**\n",
    "\n",
    "A convenient by-product of the sequential nature of the Sequential Monte Carlo (SMC) methods, is that the marginal likelihood of the model can be estimated as part of the SMC process. This is typically done as follows:\n",
    "- **Incremental Weights**: During the SMC process, at each tempering step where $ \\beta $ is increased, the weights of the samples are updated. These weights reflect the incremental importance of the samples given the change in $ \\beta $.\n",
    "- **Marginal Likelihood Estimation**: The marginal likelihood is estimated by multiplying these incremental weights across all tempering steps. This product gives an estimate of the marginal likelihood of the data under the model. Mathematically, it is often represented as the product of the average weights at each tempering step.\n",
    "- **Normalization**: These products are typically normalized at each step to prevent numerical underflow, a common issue in calculations involving products of many probabilities.\n",
    "This approach leverages the fact that SMC traverses a series of intermediate distributions, bridging the gap between the prior and the posterior, and the accumulated weights effectively capture the likelihood of the observed data across these distributions.\n",
    "\n",
    "**Motivation for using the said algorithm**:\n",
    "- It is dynamic and adaptive, which makes it highly flexible and applicable to a wide range of models, including those with complex posterior distributions that may be multi-modal or have other challenging features.\n",
    "- It incorporates the marginal likelihood calculation.\n",
    "- It is highly efficient. SMC in general is well-suited for parallel computation, making it efficient for large-scale and high-dimensional problems. This is in contrast to other high-dimensionality samplers, like Gibbs [13], where each parameter's update depends on the latest values of the other parameters, making it challenging to parallelize directly."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Results\n",
    "\n",
    "To approximate the posterior density distributions of the model parameters, we extracted the last 5000 samples of each of the final stage chains. Note that due to the dynamic nature of the used SMC algorithm, the chains grow to different lengths during the fitting process. The following figures depict the calculated posterior densities, as well as the trace paths, of all the rate and change-point location parameters for each model. Furthermore, in order to estimate the most likely values of these parameters, we computed the maximum a posteriori (MAP) estimate of their joint posterior distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "directory_path = 'plots'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model with no change-points</ins>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model0_fig2.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' not in filename and 'model0' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model0_fig3_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' in filename and 'model0' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model with one change-point</ins>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model1_fig4.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' not in filename and 'model1' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model1_fig5_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model1_fig6_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model1_fig7_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' in filename and 'model1' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model with two change-points</ins>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig8_Ts.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig9_lambdas.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' not in filename and 'model2' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig10_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig11_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig12_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig13_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model2_fig14_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' in filename and 'model2' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model with three change-points</ins>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig15_Ts.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig16_lambdas.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' not in filename and 'model3' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig17_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig18_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig19_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig20_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig21_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig22_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model3_fig23_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' in filename and 'model3' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model with four change-points</ins>**"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig24_Ts.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig25_lambdas.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' not in filename and 'model4' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig26_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig27_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig28_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig29_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig30_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig31_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig32_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig33_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots\\model4_fig34_{trace_var}_chains.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "for img_path in [\n",
    "    os.path.join(directory_path, filename) \n",
    "    for filename in os.listdir(directory_path) \n",
    "    if 'chains' in filename and 'model4' in filename\n",
    "]:\n",
    "    display(HTML(f\"<div style='text-align: center;'><img src='{img_path}'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Predictions</ins>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the data from the samples, we were able to derive predictions of the cumulative number of mining disasters. The median and 95% intervals of the predicted values were also computed. These results are portrayed in Figure 35. \n",
    "\n",
    "The plots give us an idea of how well the models with different number of change-points explain the observed data. It is deduced that the models with three and four change-points exhibit the best characteristics."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots/predictions.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(f\"<div style='text-align: center;'><img src='plots/predictions.png'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Model Selection</ins>**\n",
    "\n",
    "Occam’s razor, which is itself a consequence of Bayesian inference [14], suggests to favour simpler models over more complex ones. In our case, this means that models with more change-points, that tend to be more flexible, will be penalized by Bayesian inference. Maximization of the model evidence will indicate the balance between flexibility (which tends to improve the model's accuracy on the known data), and simplicity (which requires the model to be general enough to perform accurately in unknown cases).\n",
    "\n",
    "Figure 36 shows that the marginal likelihood considerably improves while moving from zero to three change-points, but only slightly increases when moving from three to four change-points. Thereafter, the values are expected to decrease monotonically for a larger number of change-points."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots/marginals.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(f\"<div style='text-align: center;'><img src='plots/marginals.png'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "**<ins>Convergence evaluation</ins>**"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Trace plots for the chains of all parameter paths are displayed above, in Figure 3 for the case of no change-points, Figures 5-7 for the one change-point case, Figures 10-14 for two change-points, Figures 17-23 for three change-points, and Figures 26-34 for four change-points. For each parameter, the chains appear to be mixing and converging to the same stationary distribution (i.e. the target distribution). There are no signs of issues like trend, blocking, poor mixing, or low acceptance rate.\n",
    "\n",
    "For a more objective assessment, we verified this conclusion using the Gelman-Rubin statistic $ \\hat{R} $ [15], one of the most frequently used and generally accepted approaches to monitoring convergence of MCMC samplers. An $ \\hat{R} $ value close to 1 indicates that the chains have likely converged to the same distribution (the variance within the chains is the same as the variance across the chains). Commonly, a threshold of 1.1 or 1.05 is used, where values below this threshold suggest sufficient convergence. As we can see in Figure 37, the $ \\hat{R} $ statistic remained much lower than 1.05 for all models."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div style='text-align: center;'><img src='plots/rhat.png'/></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "display(HTML(f\"<div style='text-align: center;'><img src='plots/rhat.png'/></div>\"))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Conclusion\n",
    "\n",
    "In this study, we have attempted to solve the well-known UK mining disasters changepoint detection problem, following a bayesian inference approach. We have started by presenting the problem and its dataset [5], and by formulating it as a statistical model. Next, we derived an analytical formulation to obtain the posterior distribution of the model's parameters given a variable number of change-points, which highlighted the need to employ a numerical method to approximately compute the relevant terms. To this purpose, we used a Sequential Monte Carlo sampler, which has specifically been designed to tackle high-dimensional problems efficiently [11]. The chosen algorithm, namely cascading adaptive transitional Metropolis in parallel (CATMIP), combines the Metropolis algorithm with elements of sequential importance sampling and genetic algorithms to dynamically optimize the algorithm’s efficiency as it runs [11]. We present the results of running the algorithm on the model cases of 0 - 4 change-points, using 8 parallel chains per model, and having fine-tuned one of its hyper-parameters. Furthermore, we computed the marginal likelihood of each model case, and observed that the choice of 3 or 4 change-points is optimal. Lastly, we evaluated the trace paths, as well as the $ \\hat{R} $ statistic [15], for all the sampled parameters, confirming that the used algorithm has converged in all cases."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Suggestions for future work\n",
    "\n",
    "In the process of conducting this study, it became apparent that there is a wide variety of MCMC methods, each one with its own proposed range of applications. It would certainly be beneficial to carry out a detailed literature review of all available algorithms, which was not possible in the course of this study due to time constraints. We would then experiment with multiple potential methods and compare their results to evaluate which one performs better. For a start, we would try tackling the problem with different SMC kernels, or with the Gibbs sampling approach, which, even though less efficient, might provide better convergence characteristics.\n",
    "\n",
    "Furthermore, it would be of interest to try adding an extra hierarchical level to improve our problem's statistical modeling. Specifically, it might be beneficial to sample the parameters $\\alpha$ and $\\beta$ of our Gamma priors from another prior distribution, e.g. a normal distribution. Essentially, this would drive our Bayesian framework to a higher-dimensional space, where the optimal solution might approximate the real time-series more accurately. However, this would increase the computational complexity even more, so the number of change-points we could simulate would probably be even more limited, given the available computational resources."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### References\n",
    "\n",
    "[1]\tG. Udny Yule. *VII. On a method of investigating periodicities disturbed series, with special reference to Wolfer's sunspot numbers*. Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, 226(636-646):267-298, 1927.\n",
    "\n",
    "[2]\tMontanez GD, Amizadeh S, Laptev N. *Inertial Hidden Markov Models: Modeling Change in Multivariate Time Series*. AAAI Conference on Artificial Intelligence. 2015:1819-1825.\n",
    "\n",
    "[3]\tKawahara Y, Sugiyama M. *Sequential Change-Point Detection Based on Direct Density-Ratio Estimation*. SIAM International Conference on Data Mining. 2009:389-400\n",
    "\n",
    "[4]\tPhilip Bretz. *Notes on Bayesian Changepoint Detection*. November 19, 2020.\n",
    "\n",
    "[5]\tJarrett, R.G. (1979). *A note on the intervals between coal-mining disasters*. Biometrika, 66, 191-193.\n",
    "\n",
    "[6]\tAldcroft, Derek H., *The Development of British Industry and Foreign Competition 1875-1914*. University of Toronto Press, 1968.\n",
    "\n",
    "[7]\tSalvatier, John, Thomas V. Wiecki, and Christopher Fonnesbeck. *Probabilistic programming in Python using PyMC3*. PeerJ Computer Science 2 (2016):\n",
    "\n",
    "[8]\tChristian A. Naesseth, Fredrik Lindsten and Thomas B. Schon. *Elements of Sequential Monte Carlo*. Foundations and Trends in Machine Learning, 12(3): 307-392, 2019.\n",
    "\n",
    "[9]\tDel Moral, Pierre, Arnaud Doucet, and Ajay Jasra. 2006. *Sequential Monte Carlo Samplers*. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 68 (3): 411-36.\n",
    "\n",
    "[10] Ching, J., and Chen, Y. (2007). *Transitional Markov Chain Monte Carlo Method for Bayesian Model Updating, Model Class Selection, and Model Averaging*. J. Eng. Mech., 2007, 133(7), pp. 816-832.\n",
    "\n",
    "[11] Minson, S. E., Simons, M., and Beck, J. L. (2013). *Bayesian inversion for finite fault earthquake source models I- Theory and algorithm*. Geophysical Journal International, 2013, 194(3), pp.1701 -1726.\n",
    "\n",
    "[12] Metropolis, Nicholas, et al. *Equation of State Calculations by Fast Computing Machines*. Journal of Chemical Physics, vol. 21, no. 6, 1953, pp. 1087-\n",
    "\n",
    "[13] Geman, Stuart, and Donald Geman. *Stochastic Relaxation, Gibbs Distributions, and the Bayesian Restoration of Images*. IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 6, no. 6, 1984, pp. 721-741.\n",
    "\n",
    "[14] Jefferys, William H., and James O. Berger. *Ockham’s Razor and Bayesian Analysis*. American Scientist, vol. 80, no. 1, 1992, pp. 64–72.\n",
    "\n",
    "[15] Gelman, Andrew, and Donald B. Rubin. *Inference from Iterative Simulation Using Multiple Sequences*. Statistical Science, vol. 7, no. 4, 1992, pp. 457-511."
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "pymc_env",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.11.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}