GitHub Actions build openturns/otbenchmark 10960775733

otbenchmark/master/_downloads/4314b698e56f611a280931fc186637f1/plot_check_reliability_reference_probabilities.ipynb

@@ -0,0 +1,185 @@
+{
+  "cells": [
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "\n# Check reliability problems reference probabilities\n"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "In this example, we check that the reference probabilities in the reliability problems\nare consistent with confidence bounds from Monte-Carlo simulations.\nThese 95% confidence bounds are stored in 'reliability_compute_reference_proba.csv'\nand required approximately than :math`10^9` function evaluations for each problem.\n\nWe consider two different metrics:\n\n* we check if the reference probability is within the 95% confidence bounds,\n* we compute the number of significant digits by comparing the Monte-Carlo estimator and the reference value.\n\nThe number of significant digits may be as high as 17 when all decimal digits are correct.\nHowever, the reference probabilities are only known up to 3 digits for most problems.\nIn order to keep some safeguard, we will be happy with 2 correct digits.\n\nThese metrics may fail.\n\n* On average, a fraction equal to 5% of the estimates are outside the confidence bounds.\n* The Monte-Carlo estimator may not be accurate enough, e.g. if the probability is very close to zero.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "import otbenchmark as otb\nimport pandas as pd"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "df = pd.read_csv(\"reliability_compute_reference_proba.csv\")\ndf"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "data = df.values"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "pf_reference = data[:, 1]\npmin = data[:, 3]\npmax = data[:, 4]"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "benchmarkProblemList = otb.ReliabilityBenchmarkProblemList()\nnumberOfProblems = len(benchmarkProblemList)\nnumberOfProblems"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "digitsMinimum = 2"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "categoryA = []\ncategoryB = []\ncategoryC = []\ncategoryD = []"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "for i in range(numberOfProblems):\n    problem = benchmarkProblemList[i]\n    name = problem.getName()\n    pf = problem.getProbability()\n    event = problem.getEvent()\n    antecedent = event.getAntecedent()\n    distribution = antecedent.getDistribution()\n    dimension = distribution.getDimension()\n    if pf > pmin[i] and pf < pmax[i]:\n        tagBounds = \"In\"\n    else:\n        tagBounds = \"Out\"\n    digits = otb.ComputeLogRelativeError(pf_reference[i], pf)\n    if tagBounds == \"In\" and digits >= digitsMinimum:\n        categoryA.append(name)\n    elif tagBounds == \"Out\" and digits >= digitsMinimum:\n        categoryB.append(name)\n    elif tagBounds == \"In\" and digits < digitsMinimum:\n        categoryC.append(name)\n    else:\n        categoryD.append(name)\n    print(\n        \"#%d, %-10s, pf=%.2e, ref=%.2e, C.I.=[%.2e,%.2e], digits=%d : %s\"\n        % (i, name[0:10], pf, pf_reference[i], pmin[i], pmax[i], digits, tagBounds)\n    )"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "There are four different cases.\n\n* Category A: all good. For some problems, both metrics are correct in the sense\n  that the reference probability is within the bounds and the number of significant digits is larger than 2.\n  The RP24, RP55, RP110, RP63, R-S, Axial stressed beam problems fall in that category.\n* Category B: correct digits, not in bounds.\n  We see that the RP8 problem has a reference probability outside of the 95% confidence bounds,\n  but has 2 significant digits.\n* Category C: insufficient digits, in bounds. The difficult RP28 problem fall in that category.\n* Category D: insufficient digits, not in bounds. These are suspicious problems.\n\n"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(categoryA)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(categoryB)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(categoryC)"
+      ]
+    },
+    {
+      "cell_type": "code",
+      "execution_count": null,
+      "metadata": {
+        "collapsed": false
+      },
+      "outputs": [],
+      "source": [
+        "print(categoryD)"
+      ]
+    },
+    {
+      "cell_type": "markdown",
+      "metadata": {},
+      "source": [
+        "The number of suspicious problems seems very large.\nHowever, we notice that all these cases are so that the reference probability is close,\nin absolute value, to the Monte-Carlo estimator.\n\n"
+      ]
+    }
+  ],
+  "metadata": {
+    "kernelspec": {
+      "display_name": "Python 3",
+      "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.9.20"
+    }
+  },
+  "nbformat": 4,
+  "nbformat_minor": 0
+}

otbenchmark/master/_downloads/69a98ff6834400342976e79bf25a2f75/plot_lognormal_bbrc.py

@@ -0,0 +1,156 @@
+"""
+Analysis of the LogNormal distribution in the BBRC
+==================================================
+"""
+
+# %%
+# The goal of this document is to clarify the parametrization of the LogNormal distribution in the BBRC.
+#
+# From the RPREPO
+# ---------------
+#
+# https://rprepo.readthedocs.io/en/latest/distributions.html#sec-lognormal
+#
+# * Type : univariate, continuous
+# * Support : :math:`x\in(0,\infty)`
+# * Parameter : :math:`\theta_1=\mu \in (-\infty,\infty)`, shape
+# * Parameter : :math:`\theta_2=\sigma \in (0,\infty)`, scale
+# * Mean : :math:`e^{\mu+\frac{\sigma^2}{2}}`
+# * Variance : :math:`\left(e^{\sigma^2} - 1\right) e^{2\mu + \sigma^2}`
+#
+# .. math::
+#     f(x) = \frac{1}{x} \frac{1}{\sigma \sqrt{2\pi}} \exp\left(-\frac{\left(\ln(x) - \mu\right)^2}{2 \sigma^2}\right)
+#
+# From this description we see that:
+#
+# * :math:`\theta_1` is the mean of the underlying gaussian and :math:`\theta_2`
+#   is the standard deviation of the underlying gaussian
+# * Mean is the Mean of the LogNormal random variable and Std
+#   is the standard deviation of the LogNormal random variable.
+#
+# From OpenTURNS
+# --------------
+#
+# * http://openturns.github.io/openturns/master/user_manual/_generated/openturns.LogNormal.html
+# * http://openturns.github.io/openturns/master/user_manual/_generated/openturns.LogNormalMuSigma.html
+#
+# One of the two following parametrizations must be chosen:
+#
+# * `LogNormal` with :math:`\mu_\ell=\theta_1, \sigma_\ell=\theta_2` where :math:`\theta_1` is the mean
+#    of the underlying gaussian and :math:`\theta_2` is the standard deviation of the underlying gaussian
+# * `LogNormalMuSigma` with Mean, Std where Mean is the Mean of the LogNormal random variable
+#    and Std is the standard deviation of the LogNormal random variable.
+#
+# Problem
+# -------
+#
+# The problem is to select the parametrization that best corresponds to the problem and the data.
+# The goal of this document is to make this selection clearer.
+
+# %%
+import openturns as ot
+ot.__version__
+
+# %%
+# RP60
+# ----
+#
+# https://rprepo.readthedocs.io/en/latest/reliability_problems.html#rp60
+
+# %%
+# RP60 with LogNormalMuSigma
+# --------------------------
+
+# %%
+Mean = 2200.0
+Std = 220.0
+parameters = ot.LogNormalMuSigma(Mean, Std)
+
+# %%
+X = ot.ParametrizedDistribution(parameters)
+
+# %%
+X.getMean()
+
+# %%
+X.getStandardDeviation()
+
+# %%
+# RP60 with LogNormal
+# -------------------
+
+# %%
+theta1 = 7.691
+theta2 = 0.09975
+X = ot.LogNormal(7.691, 0.09975, 0.0)
+X
+
+# %%
+X.getMean()
+
+# %%
+X.getStandardDeviation()
+
+# %%
+# RP8
+# ---
+#
+# https://rprepo.readthedocs.io/en/latest/reliability_problems.html#rp8
+
+# %%
+# RP8 with LogNormalMuSigma
+# -------------------------
+
+# %%
+Mean = 120.0
+Std = 12.0
+parameters = ot.LogNormalMuSigma(Mean, Std)
+parameters
+
+
+# %%
+X = ot.ParametrizedDistribution(parameters)
+X
+
+# %%
+X.getMean()
+
+# %%
+X.getStandardDeviation()
+
+# %%
+# RP8 with LogNormal
+# ------------------
+
+# %%
+theta1 = 4.783
+theta2 = 0.09975
+X = ot.LogNormal(4.783, 0.09975, 0.0)
+X
+
+# %%
+X.getMean()
+
+# %%
+X.getStandardDeviation()
+
+# %%
+# Conclusion
+# ----------
+#
+# We see that in the RP60 and RP8 problems, the Mean and Std parameters are exact
+# while :math:`\theta_1` and :math:`\theta_2` are given with 4 significant digits.
+# This leads to an approximation if the :math:`\theta_1` and :math:`\theta_2` parameters are used.
+#
+# This is why we choose the Mean and Std parameters as the parametrization for the BBRC.
+# This corresponds to the following code and comments:
+
+# %%
+Mean = 120.0
+Std = 12.0
+parameters = ot.LogNormalMuSigma(Mean, Std)
+X = ot.ParametrizedDistribution(parameters)
+
+# %%
+# * Mean is the Mean of the LogNormal random variable
+# * Std is the standard deviation of the LogNormal random variable.