03sim/x3p48_MC_importance.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "e353b2c8",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[<matplotlib.lines.Line2D at 0x7f82332accd0>]"
      ]
     },
     "execution_count": 1,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAXQAAAD5CAYAAAA3Os7hAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMywgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/MnkTPAAAACXBIWXMAAAsTAAALEwEAmpwYAAAjwElEQVR4nO3deXib5Znv8e9ted8dR04c29l3ssckQEtCoUACtKFMuQ5Ly3JgKDMsnc5MW2au0562055phzkz7QxbU6Atp0BKKWVNSSlryu6EJJDdCZA4q5OQzYkXSc/5Q3JiHCeRHdmPJf8+16VL1qs30i1Efrn9vM/7vOacQ0REkl+a7wJERCQxFOgiIilCgS4ikiIU6CIiKUKBLiKSIhToIiIpIv1kO5jZg8AlwE7n3IQOnjfgZ8BFwCHgOufc0pO9bv/+/d3QoUM7XbCISF+2ZMmSXc65YEfPnTTQgV8BdwEPHef5ucCo2G0mcG/s/oSGDh1KTU1NHG8vIiKtzOzj4z130iEX59xrwJ4T7DIPeMhFvQUUm1l558sUEZFTkYgx9Apgc5vHdbFtIiLSgxIR6NbBtg7XEzCzm8ysxsxq6uvrE/DWIiLSKhGBXgdUtXlcCWztaEfn3HznXLVzrjoY7HBMX0REuigRgf40cI1FnQHsc85tS8DriohIJ8QzbfFR4Bygv5nVAf8byABwzt0HLCQ6ZbGW6LTF67urWBEROb6TBrpz7sqTPO+AWxJWkYiIdInOFO0Fnlq2hU8amn2XISJJToHu2b7DLXx9wTLmL97ouxQRSXIKdM+aQxEA3qjd5bkSEUl2CnTPwpHolP0VW/ax71CL52pEJJkp0D0LRaIdunPw5sbdnqsRkWSmQPestUMHeF3DLiJyChTonrUGepop0EXk1CjQPWsN9ClVxWzc1cDWvYc9VyQiyUqB7lkoFuizR5cB6tJFpOsU6J61dujjBxXSPz9TgS4iXaZA96y1Q08PGGeN6M/rG3YTXU1BRKRzFOiehWPTFtPTjM+O7E/9gSbW7zzouSoRSUYKdM9C4Wg3HkgzPjOqPwCL12vYRUQ6T4HuWesYenpaGhXFOQwP5rF4va7mJCKdp0D3rHUMPZAWvZLfrFFB3tq4m8aWsM+yRCQJKdA9O9qhRwN99pggjS0R3vlwj8+yRCQJKdA9a9+hnzGslMz0NF5dp2EXEekcBbpnR2a5BKKBnpMZYOawfgp0Eek0BbpnRzp0syPbZo8OUrvzIHWfHPJVlogkIQW6Z+F2Qy4A54wJAvDaOk1fFJH4KdA9a52Hnp529KsYEcynojiH1zTsIiKdoED37EiHHjjaoZsZs0YHeb12Fy3hiK/SRCTJKNA9C7Wbtthq9uj+HGgK8d6mvR6qEpFkpED3rHWWS6BdoJ81sj+BNOOVtTt9lCUiSUiB7tnxOvTC7Ayqh5Tw0hoFuojER4HuWUezXFp9ftwA1mw/oOmLIhIXBbpnocixs1xanTsuehWjl9Wli0gcFOienahDHxHMZ1j/PP68WoEuIienQPfs6Dz0YwMd4NyxZby5YTcNTaGeLEtEkpAC3bNwJIIZpB0n0M8bW0ZzOKJrjYrISSnQPQtF3HG7c4DTh/WjICudFzXsIiInoUD3LBxxHY6ft8oIpDFrTJCX1u4kEtHFo0Xk+OIKdDObY2ZrzazWzO7o4PkiM3vGzJab2Uozuz7xpaamaId+4q/h8+PKqD/QxPtb9vVQVSKSjE4a6GYWAO4G5gLjgSvNbHy73W4BVjnnJgPnAP/XzDITXGtKOlmHDjB7dBlpBi+u3tFDVYlIMoqnQ58B1DrnNjrnmoEFwLx2+zigwMwMyAf2AJqWEYdQJHLCMXSAfnmZTB9Swp9WKdBF5PjiCfQKYHObx3WxbW3dBYwDtgLvA193zmmZwDjE06EDXHjaQNZsP8Cm3TprVEQ6Fk+gd5Q27Y/OXQgsAwYBU4C7zKzwmBcyu8nMasyspr5ea31DdB76yTp0iAY6wKKV27u7JBFJUvEEeh1Q1eZxJdFOvK3rgSdcVC3wITC2/Qs55+Y756qdc9XBYLCrNaeUcMR9ai3046nql8v48kIFuogcVzyB/i4wysyGxQ50XgE83W6fTcB5AGY2ABgDbExkoakqnlkurS48bSBLNn3CzgON3VyViCSjkyaJcy4E3AosAlYDjznnVprZzWZ2c2y3fwHOMrP3gReBbzvndGpjHOIdQwe4cMIAnIMXdHBURDqQHs9OzrmFwMJ22+5r8/NW4ILEltY3xDPLpdWYAQUMKc1l0codXD1zSDdXJiLJRmeKetaZDt3MmHPaQN7csIv9jS3dXJmIJBsFumcnW8ulvQtOG0hL2GmNdBE5hgLds8506ABTq4opK8ji+Q8020VEPk2B7ll0Hnr8X0NamjFnwkBeXrtTa6SLyKco0D3rbIcOcPHEchpbIryoYRcRaUOB7lkoEiE9jhOL2jp9aD/KCrJ4bkX787tEpC9ToHvWlQ49Lc24aGI5L6+t56CGXUQkRoHuWWdnubS6ZFI5zaEIf9ZJRiISo0D3rCsdOsC0wSWUF2Xz7Ipt3VCViCQjBbpnnVnLpa3WYZfX1tXrJCMRARTo3nW1Qwe4eFI5zeEIL6zUsIuIKNC968xaLu1NrSqmojiHZzXbRURQoHsXDne9QzczLplczuL1u9jT0JzgykQk2SjQPQtFXKfnobd16ZQKQhGnOekiokD37VTG0AHGlRcydmABf3hvSwKrEpFkpED3rKuzXNq6dGoFSzft5ePdDQmqSkSSkQLds1Pt0AG+OHkQZvDUMg27iPRlCnTPwl08U7StQcU5zBzWjyff24JzLkGViUiyUaB7logOHeBLUyvYuKuBFXX7ElCViCQjBbpnpzIPva05E8rJDKTx5DIdHBXpqxToHkUijoiDwCkeFAUoysngvHFlPL1sKy3hSAKqE5Fko0D3KBwb7z6VeehtfXl6JbsbmnW9UZE+SoHuUTgSDfREjKEDzB4dJFiQxe+W1CXk9UQkuSjQPQrFAj0RY+gA6YE0LptawUtrdlJ/oCkhrykiyUOB7lE4nNgOHeDy6krCEceTOnNUpM9RoHsUikQPXiYy0EeWFTB1cDGP1WzWnHSRPkaB7lGix9BbXT69ivU7D7Jcc9JF+hQFukeJHkNvdcnkcrIz0nisZnNCX1dEejcFukdHO/TEfg2F2RnMnVDOM8u2cqg5lNDXFpHeS4HuUXd16ABXzhjMgaYQzy7XRaRF+goFukfhbjgo2ur0oSWMLMvn4Xc2Jfy1RaR3UqB71J0duplx9czBLN+8lw+26OCoSF8QV6Cb2RwzW2tmtWZ2x3H2OcfMlpnZSjN7NbFlpqZQN8xDb+uyqZVkpafxiLp0kT7hpIFuZgHgbmAuMB640szGt9unGLgH+KJz7jTg8sSXmnpaD4omai2X9opyM/jC5EE89d4WDjbp4KhIqounQ58B1DrnNjrnmoEFwLx2+1wFPOGc2wTgnNPqUHEIddMsl7aumjmYhuYwT2lZXZGUF0+SVABtJzTXxba1NRooMbNXzGyJmV2TqAJTWbgbx9BbTa0qZlx5IQ+/tUlnjoqkuHgCvaO0aZ8M6cB04GLgQuA7Zjb6mBcyu8nMasyspr6+vtPFppruOPW/vdaDo6u27Wfppr3d9j4i4l88gV4HVLV5XAm0vxpxHfC8c67BObcLeA2Y3P6FnHPznXPVzrnqYDDY1ZpTRk906ACXTaugIDudX73xUbe+j4j4FU+gvwuMMrNhZpYJXAE83W6fp4CzzSzdzHKBmcDqxJaaekLdtJZLe7mZ6VxxehV/fH8b2/c1dut7iYg/Jw1051wIuBVYRDSkH3POrTSzm83s5tg+q4HngRXAO8D9zrkPuq/s1NC6fG56Nx4UbXXNmUMJO8fDb3/c7e8lIn6kx7OTc24hsLDdtvvaPb4TuDNxpaW+nurQAar65fL5cQN45O1N3PK5kWRnBLr9PUWkZ+lMUY+6ex56e9efNZTdDc08u0Lru4ikIgW6Rz0xy6WtM0eUMnpAPr98/UNNYRRJQQp0j3pqlksrM+O6s4axcut+3v5wT4+8p4j0HAW6Rz05ht7qsmkVlOZlMv+1jT32niLSMxToHh3t0Hvua8jOCHDNmUN5ac1O1u840GPvKyLdT4HukY8OHeCrZw4hOyONXyxWly6SShToHoXD0YOiPTWG3qpfXiaXT6/iyfe2snO/TjQSSRUKdI+OdOg9NG2xrRvPHkZLJMIvtRyASMpQoHvU07Nc2hpSmsec0wbym7c+1lrpIilCge6RrzH0Vl+bPYIDjSF+85aWAxBJBQp0j3zMcmlrSlUxZ4/qz/2LN3K4OeylBhFJHAW6R60duqcGHYBbPzeSXQebWfCurjsqkuwU6B6FIxHS0wwzf4k+c3gpM4b14+evbqQppC5dJJkp0D0KRZy38fO2bjt3JNv3N/L7JbruqEgyU6B7FA47LzNc2vvsyP5MrirmnldqaYnNjReR5KNA96i3dOhmxm2fG0ndJ4f5w1J16SLJSoHuUTjiSA/0jq/gvHFlTKos4mcvrqc5pC5dJBn1jjTpo3pLhw7RLv0fLhjDlr2H+W3NZt/liEgXKNA9ap3l0lvMGtWf04eWcNdL62ls0YwXkWSjQPeoN3XocLRL37G/SWePiiQhBbpH4UjvmOXS1hnDS/nsyP7c+8oGGrTGi0hSUaB71Ns69FZ/f8Fodjc088BfPvRdioh0ggLdo+g89N73FUwbXMKc0wby81c3UH+gyXc5IhKn3pcmfUhv7dABvjVnDI2hCP/14nrfpYhInBToHoUjEdI9XNwiHsOD+Vw1YzCPvLOJDfUHfZcjInFQoHvUmzt0gNvPG0V2ehp3Pr/WdykiEgcFuke9cZZLW8GCLL42ewTPr9xOzUd7fJcjIiehQPeot3foEL326IDCLH7w7CoisfXbRaR3UqB7FO3Qe/dXkJuZzh1zx7Kibh+PL6nzXY6InEDvTpMUF06CDh3g0ikVTBtczL8tWsP+xhbf5YjIcSjQPertY+itzIzvffE0djc089+axijSa8UV6GY2x8zWmlmtmd1xgv1ON7OwmX05cSWmrmQYQ281qbKYy6dX8svXP9I0RpFe6qSBbmYB4G5gLjAeuNLMxh9nv58AixJdZKrqzfPQO/LNC8eSkxHge0+vxDkdIBXpbeLp0GcAtc65jc65ZmABMK+D/W4Dfg/sTGB9KS0UcaR5vEB0ZwULsvjmnDEsXr+LZ1Zs812OiLQTT6BXAG2veFAX23aEmVUAXwLuS1xpqS9ZxtDbunrmECZVFvGDZ1ax77AOkIr0JvEEekeJ0/737Z8C33bOnfCqCGZ2k5nVmFlNfX19nCWmrlDYEejl0xbbC6QZP7p0Insamvj3RTqDVKQ3iSdN6oCqNo8rga3t9qkGFpjZR8CXgXvM7NL2L+Scm++cq3bOVQeDwa5VnEKSsUMHmFhZxDVnDuU3b3/Mss17fZcjIjHxBPq7wCgzG2ZmmcAVwNNtd3DODXPODXXODQUeB/7WOfdkootNNaGII5BEB0Xb+ocLRlNWkMUdv1+hi0qL9BInDXTnXAi4lejsldXAY865lWZ2s5nd3N0FprLedk3RzijIzuCHl05kzfYD3PNKre9yRARIj2cn59xCYGG7bR0eAHXOXXfqZfUNyTQPvSPnjx/AvCmDuOulWi48bSDjygt9lyTSpyXXEbkUk6xj6G197wunUZybwbceX0EorKEXEZ8U6B5FO/Tk/gpK8jL5wbwJvL9lHz9/baPvckT6tOROkySXCh06wEUTy7l4Yjk//fM6Vm7d57sckT5Lge6Jcy5pVluMxw8vnUBJbibf+O0yGltOeDqCiHQTBbon4djFIlKhQ4fo0Mudl09m3Y6D3KkTjkS8UKB7EooFerLOQ+/I7NFBrj1zCA/85UNer93luxyRPkeB7kmqdeit7pg7jhHBPL7x22XsPtjkuxyRPkWB7smRDj3JZ7m0l5MZ4K6rprH3cAt//9hyXYdUpAelVpokkVTt0AHGlRfy3UvG8+q6en6xWFMZRXqKAt2TUCR6Ek6qzHJp7+qZg7lo4kDuXLSWJR9/4rsckT5Bge5JKnfoEL0O6b9eNony4mxue2SpxtNFeoAC3ZNQuHUMPTUDHaAoJ4N7r57OroZmbnv0PS0NINLNFOieHOnQU2jaYkcmVBTxo0sn8MaG3dz5J81PF+lOCnRPUnWWS0cur67i6pmD+fmrG3lO1yIV6Tapnya9VKqPobf33S+MZ+rgYv7xd8u13otIN1Gge5Lqs1zay0oP8POvTKcoJ4Mbf13DzgONvksSSTkKdE/6WocOUFaYzf3XVrP3UAs3PbREi3iJJJgC3ZOjY+h9J9AhepD0P//HFJZt3ss3H1+hM0lFEkiB7snRDr3vfQVzJgzkW3PG8MzyrfybVmYUSZi4rikqidcX5qGfyN/MHsGWTw5z36sbKC/K5tqzhvouSSTpKdA96Svz0I/HzPjBvAnsPNDE955ZyYDCLOZMKPddlkhS63u/7/cSfW2WS0cCacZ/XTGVKVXF3L5gGW9oDXWRU6JA96QvznLpSE5mgAevPZ1hpXnc+FANSzdpIS+RrlKge9JXZ7l0pCQvk/93wwyCBVlc9+A7rN6233dJIklJge5JX57l0pGywmx+c8NM8rLS+eoDb7N+xwHfJYkkHaWJJ+rQj1XVL5ff3DgTM+PKX7ylUBfpJAW6J+HYQdG+Pobe3ohgPgtuOoO0WKivU6iLxE2B7klfn4d+IiOC+TzaGurz32LVVo2pi8RDge5JX5+HfjKtnXpmehpXzH+TJR/v8V2SSK+nQPdEY+gnNzyYz+9uPpPS/Cy+cv87LF5f77skkV5Nge6JZrnEp7Ikl8e+diZDSnO54Vc1PLN8q++SRHotpYkn6tDjFyzI4rc3ncmUqmJue/Q9fvHaRpzTKo0i7cUV6GY2x8zWmlmtmd3RwfNXm9mK2O0NM5uc+FJTi2a5dE5RbgYP3TCDiyeW86OFq/n+M6uO/JYjIlEnXZzLzALA3cD5QB3wrpk97Zxb1Wa3D4HZzrlPzGwuMB+Y2R0Fpwp16J2XnRHgv6+cysCibB74y4ds3nOIn14xhYLsDN+lifQK8XToM4Ba59xG51wzsACY13YH59wbzrnWRTjeAioTW2bqCYe1lktXpKUZ37lkPP8y7zReWVfPl+99k817DvkuS6RXiCfQK4DNbR7XxbYdzw3AH0+lqL5AHfqp+eqZQ/n19TPYtu8w8+5+nTc37PZdkoh38QR6R4nT4eClmX2OaKB/+zjP32RmNWZWU1/ft6eghSOOQJphpkDvqs+O6s+Tt3yG4twMvvLA2zpYKn1ePIFeB1S1eVwJHDN3zMwmAfcD85xzHbZLzrn5zrlq51x1MBjsSr0pIxQLdDk1w4P5PHXLZ7hg/AB+tHA1tzyylINNId9liXgRT6C/C4wys2FmlglcATzddgczGww8AXzVObcu8WWmnnAkovHzBCnIzuCeq6fxzxeN5fkPtnPJfy3mgy37fJcl0uNOGujOuRBwK7AIWA085pxbaWY3m9nNsd2+C5QC95jZMjOr6baKU4Q69MQyM26aNYIFN51JUyjCl+55nQf/8qGGYKRPieuaos65hcDCdtvua/PzjcCNiS0ttYUjTh16N5gxrB8Lbz+bbz6+nB88u4rF6+v5yV9Noqww23dpIt1OZ4p6ElaH3m1K8jL5xTXVfP+Lp/HGht1c+NPX+OP723yXJdLtFOieKNC7l5lx7VlDee72s6ksyeVvHl7K1xe8x56GZt+liXQbBbonoYjTwlw9YGRZPk/87Vn83edH8dyKbZz/H6/yzPKtGluXlKRE8UQdes/JCKTxd58fzbO3f5aKkhxue/Q9/vqhGuo+0RmmkloU6J6EdFC0x40dWMgTf3MW/3zRWF6v3c35//Ea976ygZZwxHdpIgmhQPckHImoQ/cgPZDGTbNG8Od/mM3Zo/rzk+fXMPdni3l1Xd8+c1lSgwLdk1BYQy4+VRTnMP+aah68rppQOMK1D77DDb96l431B32XJtJlCnRPwhGn64n2AueOHcCib8zin+aO5e0P93DBf77Gd578gPoDTb5LE+k0Bbon0TNF9Z+/N8hKD/C12SN4+R/P4YoZVTzyzibOufNl/vOFdexvbPFdnkjclCie6EzR3idYkMUPL53IC9+YxazRQX724nrO/snL3P1yLQ1a8EuSgALdk5AOivZaw4P53PuV6Txz62eZNriYOxet5ex/iwa7OnbpzRTonqhD7/0mVhbxy+tn8MTfnsXkyiLuXLSWz/z4Jf590VqNsUuvFNfiXJJ4oYgjO0OBngymDS7hl9fP4IMt+7j75VrufqWW+Ys3ctnUCm48exgjywp8lygCKNC9UYeefCZUFHHvV6azsf4gD/zlQx5fUseCdzcza3SQ684awjmjy0jTdyoeKdA9ic5D14hXMhoezOdHX5rI358/moff3sTDb3/M//xVDUNKc7lyxmAun15JaX6W7zKlD1KieKIOPfmV5mdx+3mj+Mu3z+W/r5zKgIJsfvzHNZzxry9y6yNLWby+nnBEi4BJz1GH7kkoEiGgE4tSQkYgjS9MHsQXJg9i/Y4DPPrOZn6/tI5nV2xjUFE2l02r5EvTKhgRzPddqqQ4Bbon6tBT06gBBXz3C+P51pwxvLh6J79bspl7XqnlrpdrmVhRxLwpg7hk0iAGFukKSpJ4CnRPdE3R1JadEeDiSeVcPKmcnfsbeXr5Vp5atpUfPreaHz63mtOHlnDRxHLmTBhIeVGO73IlRSjQPVGH3neUFWZz49nDufHs4WyoP8jCFdt47v1tfP+ZVXz/mVVMrizigtMGcv74AYwqy8dM/19I1yjQPdFaLn3TiGA+t503itvOG0XtzoMsWrmdP63awZ2L1nLnorVUluRw3tgyzhlbxhnDSsnJDPguWZKIAt0TdegysiyfkWUjueVzI9m+r5GX1uzkpTU7+G3NZn795sdkpqcxc1g/zh7Vn7NG9Gd8eaHmucsJKdA9CYW1loscNbAom6tmDuaqmYNpbAnzzod7eHVdPa+uq+f/LFwDQEluBmcML+WM4aXMHN6P0WUFCnj5FAW6J+rQ5XiyMwLMGh1k1ugg3wG272vkjQ27eL12N29t3M0fP9gOQHFuBtMHl1A9tB/VQ0uYWFFEdoaGaPoyBbonoYjTPHSJy8DYXPbLplUCsHnPId7auJt3P9pDzcef8OKanQCkpxnjyguZOriYSZXFTKosYkQwX78J9iEKdE/UoUtXVfXLpapfLpdXVwGw+2ATSzftZdnmT3hv015+v6SOh978GIDczADjywuZUFHE+EGFjC8vZGRZvjr5FKVA98A5p1kukjCl+VmcP34A548fAEAk4ti46yDLN+/j/S37WLl1H4/VbOZQcxiAQJoxvH8eYwYWMGZAAaMGFDBqQD5D+uWSHtD/k8lMge5B6/Ie6tClO6SlGSPLChhZVsBfTY8O00Qijo92N7B62wHWbN/P6m37WV63l2dXbDvy5zIDaQztn8uIYD7Dg3kM75/PsGAeQ0vzKMnN0Pz4JKBA9yAUiQBobFN6TFqaMTyYz/BgPhdPKj+yvaEpxPqdB6k9cjvA2u0HeGHVDkJtFhYrzE5nSGkeg/vlMrg0l6qSXKr65VBZksug4myy0jWE0xso0D1oXYFPHbr4lpeVzpSqYqZUFX9qe0s4wqY9h/hoVwMf7mrgo90NbNpzmJVb97Fo5fZPhb0ZBPOzGFScQ0VxDuVF2Qwsyqa8KIeBRVkMKMxmQGE2GRrO6XYKdA9a/zKoQ5feKiOQxohgfocrRIbCEXYcaGLznkNs3nOIrXsb2bL3EFv2Hmb1tv28uGYHjS2RY/5caV4mZYXZlBVkEYzd+udn0T8/k2B+FqX5WZTmZ1KSm6m/G12kQPcgHFaHLskrPZBGRawbP2N46THPO+fYe6iFbfsa2bG/ke37o/c7DzSxc38jO/Y3sW7HAeoPNH2q029lBsU5GZTkZdIvN5OSvExKcjMozs2kODeDktxMinMyKMrJoLDNfUFWep8/0SquQDezOcDPgABwv3Pux+2et9jzFwGHgOucc0sTXGvKONKh61dQSUFmFg3hvEzGDyo87n6RiGPf4RZ2NzRRf6CZ3Q1N7GloZtfBZvY0NPFJQwt7GprZvOcQK+qa+eRQC82hYzv/o+8L+VnpFGZnUJB99D4/O538rKO3vKyj23IzA+RlpZOTEb3PywyQkxkgNzM9KX9LOGmgm1kAuBs4H6gD3jWzp51zq9rsNhcYFbvNBO6N3UsHNIYuEj1Q2xr8I8vi+zOHm8PsPdzM3kMt7D3Uwv7GFvYdbmH/4Rb2N4Zi9y0caAxxoLGF7fsbOVgf4mBjiANNoRP+g9BeZnoauZkBcjMCZGcGyMkIkJ1x9D47I+3ofXqArDb3WekBstLTyMpIIzMQIDM9jaz0NDJjt/LYMYZEi6dDnwHUOuc2ApjZAmAe0DbQ5wEPOecc8JaZFZtZuXNu27EvJ5rlItI1OZkBcjJzuhyGzaEIh5pDHGwK0dAUpqE5REPs58Mt0ftDzSEON0c41BLicHM4emuJ3jeGos/vbmimqSVMY0uYxlAket8SJt4rDt48ewR3zB3bpc9wIvEEegWwuc3jOo7tvjvapwJIeKC///LvKH3tfwFw5L9dm/mxjuP93FZ8Qeri3K+znHO8kOkY8EoWvJURK0nhLtLdMmO34lN9ofTYrd2/K86BwxFx0b/nrnVbu58PB64C/AR6R0nT/t+hePbBzG4CbgIYPHhwHG99rPSCUjbnT+gwqs21fcs22z9VaLwX7e3ei/sGzMgYWAAZgW5/LxHpGRa7nezoWN7Aqm55/3gCvQ5o++6VwNYu7INzbj4wH6C6urpLKTau+lyoPrcrf1REJKXFM83iXWCUmQ0zs0zgCuDpdvs8DVxjUWcA+zR+LiLSs07aoTvnQmZ2K7CI6LTFB51zK83s5tjz9wELiU5ZrCU6bfH67itZREQ6Etc8dOfcQqKh3XbbfW1+dsAtiS1NREQ6Q2e2iIikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAhzzs9ZimZWD3zcxT/eH9iVwHJ80mfpnVLls6TK5wB9llZDnHPBjp7wFuinwsxqnHPVvutIBH2W3ilVPkuqfA7QZ4mHhlxERFKEAl1EJEUka6DP911AAumz9E6p8llS5XOAPstJJeUYuoiIHCtZO3QREWknKQPdzAJm9p6ZPeu7llNhZh+Z2ftmtszManzX01WxSw4+bmZrzGy1mZ3pu6auMLMxse+i9bbfzP7Od11dZWbfMLOVZvaBmT1qZtm+a+oKM/t67DOsTMbvw8weNLOdZvZBm239zOwFM1sfuy9JxHslZaADXwdW+y4iQT7nnJuS5NOxfgY875wbC0wmSb8b59za2HcxBZhOdCnoP/itqmvMrAK4Hah2zk0guvT1FX6r6jwzmwD8NdFrG08GLjGzUX6r6rRfAXPabbsDeNE5Nwp4Mfb4lCVdoJtZJXAxcL/vWgTMrBCYBTwA4Jxrds7t9VpUYpwHbHDOdfXkt94gHcgxs3Qglw6uIpYExgFvOecOOedCwKvAlzzX1CnOudeAPe02zwN+Hfv518CliXivpAt04KfAt4CI5zoSwQF/MrMlseutJqPhQD3wy9gw2P1mlue7qAS4AnjUdxFd5ZzbAvw7sInoxdr3Oef+5LeqLvkAmGVmpWaWS/RCOt1zQc6eNaD1qm6x+7JEvGhSBbqZXQLsdM4t8V1LgnzGOTcNmAvcYmazfBfUBenANOBe59xUoIEE/froS+xSi18Efue7lq6KjcnOA4YBg4A8M/uK36o6zzm3GvgJ8ALwPLAcCHktqhdLqkAHPgN80cw+AhYA55rZb/yW1HXOua2x+51Ex2pn+K2oS+qAOufc27HHjxMN+GQ2F1jqnNvhu5BT8HngQ+dcvXOuBXgCOMtzTV3inHvAOTfNOTeL6NDFet81JcAOMysHiN3vTMSLJlWgO+f+yTlX6ZwbSvRX4pecc0nXdQCYWZ6ZFbT+DFxA9NfLpOKc2w5sNrMxsU3nAas8lpQIV5LEwy0xm4AzzCzXzIzo95KUB6vNrCx2Pxi4jOT/bgCeBq6N/Xwt8FQiXjSua4pKtxgA/CH6d4104BHn3PN+S+qy24CHY0MVG0nii4THxmnPB77mu5ZT4Zx728weB5YSHaJ4j+Q90/L3ZlYKtAC3OOc+8V1QZ5jZo8A5QH8zqwP+N/Bj4DEzu4HoP76XJ+S9dKaoiEhqSKohFxEROT4FuohIilCgi4ikCAW6iEiKUKCLiKQIBbqISIpQoIuIpAgFuohIivj/b8IKOaPlT3sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.stats as stats\n",
    "import matplotlib.pyplot as plt\n",
    "import numpy as np\n",
    "\n",
    "x = np.linspace(4, 10, 100)\n",
    "plt.plot(x, stats.expon(5).pdf(x))\n",
    "plt.plot(x, stats.norm().pdf(x))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b6b2f8d8",
   "metadata": {},
   "source": [
    "The true probability of the extreme event, $\\mathrm{P} (X>5) $ is extremely small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "83684c29",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "2.866515719235352e-07"
      ]
     },
     "execution_count": 2,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "P_true = 1 - stats.norm().cdf(5)\n",
    "P_true"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "00d88ac7",
   "metadata": {},
   "source": [
    "First try ordinary MC."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "0417cd5c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(0.0, 1.0)"
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N= 10000\n",
    "y = stats.norm().rvs(N)\n",
    "P_mc = 1.0/N * np.sum(y > 5)\n",
    "# estimate and relative error\n",
    "P_mc, np.abs(P_mc - P_true)/P_true"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "3c935d65",
   "metadata": {},
   "source": [
    "The result is completely off. Now we apply importance sampling with the truncated exponential"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "cbe78b8d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(2.946190959301109e-07, 0.02779515197879704)"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "N = 10000\n",
    "y = stats.expon(loc=5).rvs(N)\n",
    "P_is = 1.0/N * np.sum(stats.norm().pdf(y)/stats.expon(loc=5).pdf(y))\n",
    "# estimate and relative error\n",
    "P_is, np.abs(P_is- P_true)/P_true"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "356f788d",
   "metadata": {},
   "source": [
    "The result is excellent, with a relative error of only $0.4\\%.$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "38b64f5e",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.7"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}