ex1_linear_limit_state_function.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from rips import *\n",
    "from rips.utils import FeynmanKac, Particle"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 1. Setting up model\n",
    "\n",
    "Here we consider the stochastic system $g(X) = \\beta  - \\frac{1}{\\sqrt{n}} \\sum_{i=1}^{n} X_i$, where $(X_i)$ are i.i.d. standard normal random variables.\n",
    "\n",
    "We are interested in the probability $p = P\\{g(X)\\leq 0\\}$.\n",
    "\n",
    "Note that when $\\beta \\to \\infty$ we have that $p \\to 0$, so crude Monte Carlo samling becomes unfeasible for large $\\beta$. \n",
    "\n",
    "Next we show how to use pims to approximate the probability $p$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Define model with linear limit state function\n",
    "class LinearLSF(FeynmanKac, StandardGaussian):\n",
    "    def __init__(self, beta=1, n=10, **kwargs):\n",
    "        self.beta = beta\n",
    "        self.n = n\n",
    "        super().__init__(**kwargs)\n",
    "    \n",
    "    def score_function(self, particle: Particle) -> float:\n",
    "        particle.response = particle.path.sum() / np.sqrt(self.n)\n",
    "        return particle.response / self.beta\n",
    "    \n",
    "    @property\n",
    "    def num_variables(self) -> int:\n",
    "        return self.n\n",
    "\n",
    "model = LinearLSF(beta=6, n=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sample a particle\n",
    "particle = model.sample_particle()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "-0.13792925318289992"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# This is the score function evaluated on a random vector\n",
    "particle.score  # Greater than one signifies event of interest g(X) <= 0"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 2. Probabilistic inference"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Adaptive levels...\n",
      "[-inf, 0.13891543474600263]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754, 0.9119740613454691]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754, 0.9119740613454691, 0.956230929679357]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754, 0.9119740613454691, 0.956230929679357, 0.999729193931889]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754, 0.9119740613454691, 0.956230929679357, 0.999729193931889, 1.041255440521729]\n",
      "[-inf, 0.13891543474600263, 0.2856820035665473, 0.4010426081789563, 0.49081591352523235, 0.569599038921804, 0.6408751412681312, 0.7008175938164891, 0.7571157746221547, 0.8115499537757305, 0.8630944693669754, 0.9119740613454691, 0.956230929679357, 0.999729193931889, 1]\n",
      "Fixed levels...\n",
      "done\n"
     ]
    },
    {
     "data": {
      "text/plain": [
       "{'cpu_time_smc': 12.560560375,\n",
       " 'p_smc': 8.0592896e-10,\n",
       " 'ng_smc': 140001,\n",
       " 'cpu_time_gs': 12.469817374999998,\n",
       " 'p_bar': 9.387212800000002e-10,\n",
       " 'ng_gs': 146950,\n",
       " 'k': 14,\n",
       " 'var': 1.8415964470528175e-16,\n",
       " 'std': 1.3570543272296867e-08}"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Number of samples\n",
    "model.num_of_particles = 10000\n",
    "# Markov Chain Monte Carlo methods\n",
    "kernels = [\n",
    "    PCN(),\n",
    "    # MMH()\n",
    "]\n",
    "# Sample-based Monte Carlo integration\n",
    "sample = ComboUQ(model, kernels)\n",
    "sample.summary_results()\n"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 3. Compare approximations to exact solution"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.865876450376946e-10"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "from scipy.stats import norm\n",
    "p_exact = norm.cdf(-6)\n",
    "p_exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "9.387212800000002e-10"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Unbiased estimator\n",
    "sample.p_bar"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.04851709351769308"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# relative error\n",
    "abs(sample.p_bar - p_exact) / p_exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8.0592896e-10"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# Biased estimator\n",
    "sample.p_smc"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.18311468418073704"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# relative error\n",
    "abs(sample.p_smc - p_exact) / p_exact"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "robspy",
   "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.6"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}