.. DO NOT EDIT. .. THIS FILE WAS AUTOMATICALLY GENERATED BY SPHINX-GALLERY. .. TO MAKE CHANGES, EDIT THE SOURCE PYTHON FILE: .. "auto_calibration/bayesian_calibration/plot_rwmh_python_distribution.py" .. LINE NUMBERS ARE GIVEN BELOW. .. only:: html .. note:: :class: sphx-glr-download-link-note :ref:`Go to the end ` to download the full example code .. rst-class:: sphx-glr-example-title .. _sphx_glr_auto_calibration_bayesian_calibration_plot_rwmh_python_distribution.py: Posterior sampling using a PythonDistribution ============================================= .. GENERATED FROM PYTHON SOURCE LINES 6-61 In this example we are going to show how to do Bayesian inference using the :class:`~openturns.RandomWalkMetropolisHastings` algorithm in a statistical model defined through a :class:`~openturns.PythonDistribution`. This method is illustrated on a simple lifetime study test-case, which involves censored data, as described hereafter. In the following, we assume that the lifetime :math:`T_i` of an industrial component follows the Weibull distribution :math:`\mathcal W(\alpha, \beta)`, with CDF given by :math:`F(t|\alpha,\beta)= 1 - e^{-\left( \frac{t}{\beta} \right)^\alpha}`. Our goal is to estimate the model parameters :math:`\alpha, \beta` based on a dataset of recorded failures :math:`(t_1, \ldots, t_n),` some of which correspond to actual failures, and the remaining are right-censored. Let :math:`(f_1, \ldots, f_n) \in \{0,1\}^n` represent the nature of each datum, :math:`f_i=1` if :math:`t_i` corresponds to an actual failure, :math:`f_i=0` if it is right-censored. Note that the likelihood of each recorded failure is given by the Weibull density: .. math:: \mathcal L(t_i | f_i=1, \alpha, \beta) = \frac{\alpha}{\beta}\left( \frac{t_i}{\beta} \right)^{\alpha-1} e^{-\left( \frac{t_i}{\beta} \right)^\alpha}. On the other hand, the likelihood of each right-censored observation is given by: .. math:: \mathcal L(t_i | f_i=0, \alpha, \beta) = e^{-\left( \frac{t_i}{\beta} \right)^\alpha}. Furthermore, assume that the prior information available on :math:`\alpha, \beta` is represented by independent prior laws, whose respective densities are denoted by :math:`\pi(\alpha)` and :math:`\pi(\beta).` The posterior distribution of :math:`(\alpha, \beta)` represents the update of the prior information on :math:`(\alpha, \beta)` given the dataset. Its PDF is known up to a multiplicative constant: .. math:: \pi(\alpha, \beta | (t_1, f_1), \ldots, (t_n, f_n) ) \propto \pi(\alpha)\pi(\beta) \left(\frac{\alpha}{\beta}\right)^{\sum_i f_i} \left(\prod_{f_i = 1} \frac{t_i}{\beta}\right)^{\alpha-1} \exp\left[-\sum_{i=1}^n\left(\frac{t_i}{\beta}\right)^\alpha\right]. The :class:`~openturns.RandomWalkMetropolisHastings` class can be used to sample from the posterior distribution. It relies on the following objects: - The prior probability density :math:`\pi(\vect{\theta})` reflects beliefs about the possible values of :math:`\vect{\theta} = (\alpha, \beta)` before the experimental data are considered. - Initial values :math:`\vect{\theta}_0` of the parameters. - An proposal distribution used to update parameters. Additionnaly we want to define the likelihood term defined by these objects: - The conditional density :math:`p(t_{1:n}|f_{1:n}, \alpha, \beta)` will be defined as a :class:`~openturns.PythonDistribution`. - The sample of observations acting as the parameters of the conditional density Set up the PythonDistribution ----------------------------- The censured Weibuill likelihood is outside the usual catalog of probability distributions in OpenTURNS, hence we need to define it using the :class:`~openturns.PythonDistribution` class. .. GENERATED FROM PYTHON SOURCE LINES 64-72 .. code-block:: Python import numpy as np import openturns as ot from openturns.viewer import View ot.Log.Show(ot.Log.NONE) ot.RandomGenerator.SetSeed(123) .. GENERATED FROM PYTHON SOURCE LINES 73-82 The following methods must be defined: - `getRange`: required for conversion to the :class:`~openturns.Distribution` format - `computeLogPDF`: used by :class:`~openturns.RandomWalkMetropolisHastings` to evaluate the posterior density - `setParameter` used by :class:`~openturns.RandomWalkMetropolisHastings` to test new parameter values .. note:: We formally define a bivariate distribution on the :math:`(t_i, f_i)` couple, even though :math:`f_i` has no distribution (it is simply a covariate). This is not an issue, since the sole purpose of this :class:`~openturns.PythonDistribution` object is to pass the likelihood calculation over to :class:`~openturns.RandomWalkMetropolisHastings`. .. GENERATED FROM PYTHON SOURCE LINES 84-120 .. code-block:: Python class CensoredWeibull(ot.PythonDistribution): """ Right-censored Weibull log-PDF calculation Each data point x is assumed 2D, with: x[0]: observed functioning time x[1]: nature of x[0]: if x[1]=0: x[0] is a censoring time if x[1]=1: x[0] is a time-to failure """ def __init__(self, beta=5000.0, alpha=2.0): super(CensoredWeibull, self).__init__(2) self.beta = beta self.alpha = alpha def getRange(self): return ot.Interval([0, 0], [1, 1], [True] * 2, [False, True]) def computeLogPDF(self, x): if not (self.alpha > 0.0 and self.beta > 0.0): return -np.inf log_pdf = -((x[0] / self.beta) ** self.alpha) log_pdf += (self.alpha - 1) * np.log(x[0] / self.beta) * x[1] log_pdf += np.log(self.alpha / self.beta) * x[1] return log_pdf def setParameter(self, parameter): self.beta = parameter[0] self.alpha = parameter[1] def getParameter(self): return [self.beta, self.alpha] .. GENERATED FROM PYTHON SOURCE LINES 121-122 Convert to :class:`~openturns.Distribution` .. GENERATED FROM PYTHON SOURCE LINES 124-129 .. code-block:: Python conditional = ot.Distribution(CensoredWeibull()) .. GENERATED FROM PYTHON SOURCE LINES 130-134 Observations, prior, initial point and proposal distributions ------------------------------------------------------------- Define the observations .. GENERATED FROM PYTHON SOURCE LINES 137-143 .. code-block:: Python Tobs = np.array([4380, 1791, 1611, 1291, 6132, 5694, 5296, 4818, 4818, 4380]) fail = np.array([True] * 4 + [False] * 6) x = ot.Sample(np.vstack((Tobs, fail)).T) .. GENERATED FROM PYTHON SOURCE LINES 144-146 Define a uniform prior distribution for :math:`\alpha` and a Gamma prior distribution for :math:`\beta`. .. GENERATED FROM PYTHON SOURCE LINES 148-160 .. code-block:: Python alpha_min, alpha_max = 0.5, 3.8 a_beta, b_beta = 2, 2e-4 priorCopula = ot.IndependentCopula(2) # prior independence priorMarginals = [] # prior marginals priorMarginals.append(ot.Gamma(a_beta, b_beta)) # Gamma prior for beta priorMarginals.append(ot.Uniform(alpha_min, alpha_max)) # uniform prior for alpha prior = ot.ComposedDistribution(priorMarginals, priorCopula) prior.setDescription(["beta", "alpha"]) .. GENERATED FROM PYTHON SOURCE LINES 161-163 We select prior means as the initial point of the Metropolis-Hastings algorithm. .. GENERATED FROM PYTHON SOURCE LINES 165-168 .. code-block:: Python initialState = [a_beta / b_beta, 0.5 * (alpha_max - alpha_min)] .. GENERATED FROM PYTHON SOURCE LINES 169-172 For our random walk proposal distributions, we choose normal steps, with standard deviation equal to roughly :math:`10\%` of the prior range (for the uniform prior) or standard deviation (for the normal prior). .. GENERATED FROM PYTHON SOURCE LINES 174-180 .. code-block:: Python proposal = [] proposal.append(ot.Normal(0.0, 0.1 * np.sqrt(a_beta / b_beta ** 2))) proposal.append(ot.Normal(0.0, 0.1 * (alpha_max - alpha_min))) proposal = ot.ComposedDistribution(proposal) .. GENERATED FROM PYTHON SOURCE LINES 181-183 Sample from the posterior distribution -------------------------------------- .. GENERATED FROM PYTHON SOURCE LINES 185-193 .. code-block:: Python sampler = ot.RandomWalkMetropolisHastings(prior, initialState, proposal) sampler.setLikelihood(conditional, x) sampleSize = 10000 sample = sampler.getSample(sampleSize) # compute acceptance rate print("Acceptance rate: %s" % (sampler.getAcceptanceRate())) .. rst-class:: sphx-glr-script-out .. code-block:: none Acceptance rate: 0.3347 .. GENERATED FROM PYTHON SOURCE LINES 194-196 Plot prior to posterior marginal plots. .. GENERATED FROM PYTHON SOURCE LINES 198-211 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) grid = ot.GridLayout(1, 2) grid.setTitle("Bayesian inference") for parameter_index in range(2): graph = posterior.getMarginal(parameter_index).drawPDF() priorGraph = prior.getMarginal(parameter_index).drawPDF() graph.add(priorGraph) graph.setColors(ot.Drawable.BuildDefaultPalette(2)) graph.setLegends(["Posterior", "Prior"]) grid.setGraph(0, parameter_index, graph) _ = View(grid) .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_rwmh_python_distribution_001.png :alt: Bayesian inference :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_rwmh_python_distribution_001.png :class: sphx-glr-single-img .. GENERATED FROM PYTHON SOURCE LINES 212-214 Define an improper prior -------------------------- .. GENERATED FROM PYTHON SOURCE LINES 216-221 Now, define an improper prior: .. math:: \mathcal \pi(\beta, \alpha) \propto \frac{1}{\beta}. .. GENERATED FROM PYTHON SOURCE LINES 221-226 .. code-block:: Python logpdf = ot.SymbolicFunction(["beta", "alpha"], ["-log(beta)"]) support = ot.Interval([0] * 2, [1] * 2) support.setFiniteUpperBound([False] * 2) .. GENERATED FROM PYTHON SOURCE LINES 227-228 Sample from the posterior distribution .. GENERATED FROM PYTHON SOURCE LINES 228-234 .. code-block:: Python sampler2 = ot.RandomWalkMetropolisHastings(logpdf, support, initialState, proposal) sampler2.setLikelihood(conditional, x) sample2 = sampler2.getSample(1000) print("Acceptance rate: %s" % (sampler2.getAcceptanceRate())) .. rst-class:: sphx-glr-script-out .. code-block:: none Acceptance rate: 0.37 .. GENERATED FROM PYTHON SOURCE LINES 235-237 Plot posterior marginal plots only as prior cannot be drawn meaningfully. .. GENERATED FROM PYTHON SOURCE LINES 239-251 .. code-block:: Python kernel = ot.KernelSmoothing() posterior = kernel.build(sample) grid = ot.GridLayout(1, 2) grid.setTitle("Bayesian inference (with log-pdf)") for parameter_index in range(2): graph = posterior.getMarginal(parameter_index).drawPDF() graph.setColors(ot.Drawable.BuildDefaultPalette(2)) graph.setLegends(["Posterior"]) grid.setGraph(0, parameter_index, graph) _ = View(grid) View.ShowAll() .. image-sg:: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_rwmh_python_distribution_002.png :alt: Bayesian inference (with log-pdf) :srcset: /auto_calibration/bayesian_calibration/images/sphx_glr_plot_rwmh_python_distribution_002.png :class: sphx-glr-single-img .. rst-class:: sphx-glr-timing **Total running time of the script:** (0 minutes 5.505 seconds) .. _sphx_glr_download_auto_calibration_bayesian_calibration_plot_rwmh_python_distribution.py: .. only:: html .. container:: sphx-glr-footer sphx-glr-footer-example .. container:: sphx-glr-download sphx-glr-download-jupyter :download:`Download Jupyter notebook: plot_rwmh_python_distribution.ipynb ` .. container:: sphx-glr-download sphx-glr-download-python :download:`Download Python source code: plot_rwmh_python_distribution.py `