thirteen_fitting_arch_compressive_failure_probability_data_stack_implementation.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 67,
   "metadata": {},
   "outputs": [],
   "source": [
    "#https://stackoverflow.com/questions/73071058/fit-data-with-a-lognormal-function-via-maximum-likelihood-estimators/73112103#73112103\n",
    "#https://towardsdatascience.com/log-normal-distribution-a-simple-explanation-7605864fb67c"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 68,
   "metadata": {},
   "outputs": [],
   "source": [
    "from functools import partial\n",
    "import numpy as np\n",
    "from scipy import optimize, stats\n",
    "import matplotlib.pyplot as plt"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 69,
   "metadata": {},
   "outputs": [],
   "source": [
    "# %store"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 70,
   "metadata": {},
   "outputs": [],
   "source": [
    "%store -r pier_local_scour_depth\n",
    "%store -r arch_compressive_failure_points_number\n",
    "%store -r convergence_simulation_number"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 71,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Failure probability data from Baker (2015), https://doi.org/10.1193/021113EQS025M\n",
    "# im = np.array([0.2, 0.3, 0.4, 0.6, 0.7, 0.8, 0.9, 1])\n",
    "# number_of_analyses = np.array([40, 40, 40, 40, 40, 40, 40, 40])\n",
    "# number_of_collapses = np.array([0, 0, 0, 4, 6, 13, 12, 16])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 72,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Failure probability data from Qeshta et al., 2021, https://doi.org/10.1080/15732479.2021.1892774\n",
    "# im = np.array([4.5, 5, 5.5, 6, 6.5, 7, 7.5, 8, 8.5, 9, 9.46])\n",
    "# number_of_analyses = np.array([10, 10, 10, 10, 10, 10, 10, 10, 10, 10, 10]) \n",
    "# number_of_collapses = np.array([0, 0, 0, 0, 0, 0, 1, 3, 7, 9, 10])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 89,
   "metadata": {},
   "outputs": [],
   "source": [
    "im = []\n",
    "\n",
    "#if im has an initial value of 0, this leads to -infinite when log is calculated with the lognormal pdf. Therefore, 0 is approximated with 10^-6 for which the code still works.\n",
    "\n",
    "def replace_zero_im(intensity_measure_vector):\n",
    "    \n",
    "    global im #https://scipy-lectures.org/intro/language/functions.html\n",
    "\n",
    "    for i in range(len(intensity_measure_vector)):\n",
    "        if intensity_measure_vector[i]==0:\n",
    "            intensity_measure_vector[i]=10**-6\n",
    "    \n",
    "    im = intensity_measure_vector\n",
    "    return im"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 94,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Failure probability data from \"twelve_arch_compressive_limit_state_equation.ipynb\"\n",
    "im_input = np.array([0, 0.289, 0.578, 0.867, 1.156, 1.444, 1.733, 2.022, 2.311, 2.6])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 95,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "im =  [1.000e-06 2.890e-01 5.780e-01 8.670e-01 1.156e+00 1.444e+00 1.733e+00\n",
      " 2.022e+00 2.311e+00 2.600e+00]\n"
     ]
    }
   ],
   "source": [
    "replace_zero_im(im_input)\n",
    "\n",
    "print('im = ', im) #https://stackoverflow.com/questions/21610198/runtimewarning-divide-by-zero-encountered-in-log\n",
    "number_of_analyses = np.array([100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000])\n",
    "number_of_collapses = np.array([69, 69, 69, 69, 259, 1880, 15286, 74760, 99995, 100000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 96,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Above data without the initial zero im value\n",
    "# im = np.array([0.289, 0.578, 0.867, 1.156, 1.444, 1.733, 2.022, 2.311, 2.6]) #https://stackoverflow.com/questions/21610198/runtimewarning-divide-by-zero-encountered-in-log\n",
    "# number_of_analyses = np.array([100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000, 100000])\n",
    "# number_of_collapses = np.array([69, 69, 69, 259, 1880, 15286, 74760, 99995, 100000])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 97,
   "metadata": {},
   "outputs": [],
   "source": [
    "#Failure probability data from \"twelve_arch_compressive_limit_state_equation.ipynb\" and other stored values\n",
    "# im = np.array(pier_local_scour_depth[0])\n",
    "# number_of_analyses = np.ones(len(pier_local_scour_depth[0]))*convergence_simulation_number\n",
    "# number_of_collapses = np.array(arch_compressive_failure_points_number)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 98,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "<matplotlib.collections.PathCollection at 0x26d398185e0>"
      ]
     },
     "execution_count": 98,
     "metadata": {},
     "output_type": "execute_result"
    },
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8qNh9FAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAfkUlEQVR4nO3df2zU9eHH8de1pXfq2nPAuBaptTrFsmbOXge22i3qOCmmCYmJNUZAh4tlKimdOiqJDGLSbHOMbdoiEzRG5hoFjGQNcsnkh4Jh1CMb1qkZ1Va52hSWu4qjSPv5/sG3jedd4T5Hr++74/lI7o9+eH/a9/vjJ+nTz33uU4dlWZYAAAAMyTI9AQAAcGEjRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGBUjukJxGN4eFhHjx5VXl6eHA6H6ekAAIA4WJalgYEBTZ8+XVlZY1//SIsYOXr0qIqKikxPAwAAJKCnp0czZswY89/TIkby8vIknVlMfn6+4dkAAIB4hMNhFRUVjf4eH0taxMjIWzP5+fnECAAAaeZct1hwAysAADCKGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRafHQMwAAUsnQsKUDXcfVN3BS0/Jcml0yWdlZ6fe301JlHbZjZM+ePfrtb3+rjo4OBYNBbdu2TQsWLDjrPrt371ZjY6Pee+89TZ8+XY899pjq6+sTnTMAAMbsOBzU6u2dCoZOjm4rdLu0qnaW5pUVGpyZPam0Dttv05w4cULXXXednn766bjGd3V1af78+aqurlYgENDjjz+uZcuWacuWLbYnCwCASTsOB7X0pXcjfoFLUm/opJa+9K52HA4ampk9qbYO21dGampqVFNTE/f49evX6/LLL9e6deskSaWlpTp48KCeeuop3XHHHXZ/PAAARgwNW1q9vVNWjH+zJDkkrd7eqbmzClL6LZtUXEfSb2Ddv3+/fD5fxLbbbrtNBw8e1FdffRVzn8HBQYXD4YgXAAAmHeg6HnUl4essScHQSR3oOj5xk0pAKq4j6THS29srj8cTsc3j8ej06dPq7++PuU9zc7Pcbvfoq6ioKNnTBADgrPoGxv4Fnsg4U1JxHRPy0d5v/ulgy7Jibh/R1NSkUCg0+urp6Un6HAEAOJtpea5xHWdKKq4j6R/tLSgoUG9vb8S2vr4+5eTkaMqUKTH3cTqdcjqdyZ4aAABxm10yWYVul3pDJ2Peb+GQVOA+8/HYVJaK60j6lZHKykr5/f6IbTt37lRFRYUmTZqU7B8PAMC4yM5yaFXtLElnfmF/3cjXq2pnpfTNq1JqrsN2jHzxxRc6dOiQDh06JOnMR3cPHTqk7u5uSWfeYlm0aNHo+Pr6en3yySdqbGzU+++/r02bNmnjxo165JFHxmcFAABMkHllhWq9p1wF7si3MArcLrXeU542zxlJtXU4rJEbOOK0a9cu3XzzzVHbFy9erBdeeEH33nuvPv74Y+3atWv033bv3q3ly5ePPvTsl7/8pa2HnoXDYbndboVCIeXn59uZLgAA4y5Vnlx6vpK9jnh/f9uOEROIEQAA0k+8v7/5Q3kAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMCohGKkpaVFJSUlcrlc8nq92rt371nHb968Wdddd50uvvhiFRYW6r777tOxY8cSmjAAAMgstmOkra1NDQ0NWrlypQKBgKqrq1VTU6Pu7u6Y49966y0tWrRIS5Ys0XvvvadXXnlF//jHP3T//fef9+QBAED6sx0ja9eu1ZIlS3T//fertLRU69atU1FRkVpbW2OOf+edd3TFFVdo2bJlKikp0U033aQHHnhABw8ePO/JAwCA9GcrRk6dOqWOjg75fL6I7T6fT/v27Yu5T1VVlT799FO1t7fLsix9/vnnevXVV3X77beP+XMGBwcVDocjXgAAIDPZipH+/n4NDQ3J4/FEbPd4POrt7Y25T1VVlTZv3qy6ujrl5uaqoKBAl156qf70pz+N+XOam5vldrtHX0VFRXamCQAA0khCN7A6HI6Iry3Lito2orOzU8uWLdMTTzyhjo4O7dixQ11dXaqvrx/z+zc1NSkUCo2+enp6EpkmAABIAzl2Bk+dOlXZ2dlRV0H6+vqirpaMaG5u1o033qhHH31UkvT9739fl1xyiaqrq/Xkk0+qsLAwah+n0ymn02lnagAAIE3ZujKSm5srr9crv98fsd3v96uqqirmPl9++aWysiJ/THZ2tqQzV1QAAMCFzfbbNI2NjXruuee0adMmvf/++1q+fLm6u7tH33ZpamrSokWLRsfX1tZq69atam1t1ZEjR/T2229r2bJlmj17tqZPnz5+KwEAAGnJ1ts0klRXV6djx45pzZo1CgaDKisrU3t7u4qLiyVJwWAw4pkj9957rwYGBvT000/rF7/4hS699FLdcsst+vWvfz1+qwAAAGnLYaXBeyXhcFhut1uhUEj5+fmmpwMAAOIQ7+9v/jYNAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKNsf7QXAIBEDQ1bOtB1XH0DJzUtz6XZJZOVnRX7z4ngwkGMAAAmxI7DQa3e3qlg6OTotkK3S6tqZ2leWfSfBsGFg7dpAABJt+NwUEtfejciRCSpN3RSS196VzsOBw3NDKmAGAEAJNXQsKXV2zsV6wmbI9tWb+/U0HDKP4MTSUKMAACS6kDX8agrIl9nSQqGTupA1/GJmxRSCjECAEiqvoGxQySRccg8xAgAIKmm5bnGdRwyDzECAEiq2SWTVeh2aawP8Dp05lM1s0smT+S0kEKIEQBAUmVnObSqdpYkRQXJyNeramfxvJELGDECAEi6eWWFar2nXAXuyLdiCtwutd5TznNGLnA89AwAMCHmlRVq7qwCnsCKKMQIAGDCZGc5VHnVFNPTQIrhbRoAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYlVCMtLS0qKSkRC6XS16vV3v37j3r+MHBQa1cuVLFxcVyOp266qqrtGnTpoQmDAAAMkuO3R3a2trU0NCglpYW3XjjjXr22WdVU1Ojzs5OXX755TH3ufPOO/X5559r48aN+u53v6u+vj6dPn36vCcPAADSn8OyLMvODnPmzFF5eblaW1tHt5WWlmrBggVqbm6OGr9jxw7dddddOnLkiCZPnpzQJMPhsNxut0KhkPLz8xP6HgAAYGLF+/vb1ts0p06dUkdHh3w+X8R2n8+nffv2xdzn9ddfV0VFhX7zm9/osssu0zXXXKNHHnlE//vf/8b8OYODgwqHwxEvAACQmWy9TdPf36+hoSF5PJ6I7R6PR729vTH3OXLkiN566y25XC5t27ZN/f39+vnPf67jx4+Ped9Ic3OzVq9ebWdqAAAgTSV0A6vD4Yj42rKsqG0jhoeH5XA4tHnzZs2ePVvz58/X2rVr9cILL4x5daSpqUmhUGj01dPTk8g0AQBAGrB1ZWTq1KnKzs6OugrS19cXdbVkRGFhoS677DK53e7RbaWlpbIsS59++qmuvvrqqH2cTqecTqedqQEAgDRl68pIbm6uvF6v/H5/xHa/36+qqqqY+9x44406evSovvjii9FtH374obKysjRjxowEpgwAADKJ7bdpGhsb9dxzz2nTpk16//33tXz5cnV3d6u+vl7SmbdYFi1aNDr+7rvv1pQpU3Tfffeps7NTe/bs0aOPPqqf/vSnuuiii8ZvJQAAIC3Zfs5IXV2djh07pjVr1igYDKqsrEzt7e0qLi6WJAWDQXV3d4+O/9a3viW/36+HH35YFRUVmjJliu688049+eST47cKAACQtmw/Z8QEnjMCAED6ScpzRgAAAMYbMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjEoqRlpYWlZSUyOVyyev1au/evXHt9/bbbysnJ0c/+MEPEvmxAAAgA9mOkba2NjU0NGjlypUKBAKqrq5WTU2Nuru7z7pfKBTSokWLdOuttyY8WQAAkHkclmVZdnaYM2eOysvL1draOrqttLRUCxYsUHNz85j73XXXXbr66quVnZ2t1157TYcOHYr7Z4bDYbndboVCIeXn59uZLgAAMCTe39+2roycOnVKHR0d8vl8Edt9Pp/27ds35n7PP/+8/vOf/2jVqlVx/ZzBwUGFw+GIFwAAyEy2YqS/v19DQ0PyeDwR2z0ej3p7e2Pu89FHH2nFihXavHmzcnJy4vo5zc3Ncrvdo6+ioiI70wQAAGkkoRtYHQ5HxNeWZUVtk6ShoSHdfffdWr16ta655pq4v39TU5NCodDoq6enJ5FpAgCANBDfpYr/N3XqVGVnZ0ddBenr64u6WiJJAwMDOnjwoAKBgB566CFJ0vDwsCzLUk5Ojnbu3Klbbrklaj+n0ymn02lnagAAIE3ZujKSm5srr9crv98fsd3v96uqqipqfH5+vv71r3/p0KFDo6/6+nrNnDlThw4d0pw5c85v9gAAIO3ZujIiSY2NjVq4cKEqKipUWVmpDRs2qLu7W/X19ZLOvMXy2Wef6cUXX1RWVpbKysoi9p82bZpcLlfUdgAAcGGyHSN1dXU6duyY1qxZo2AwqLKyMrW3t6u4uFiSFAwGz/nMEQAAgBG2nzNiAs8ZAQAg/STlOSMAAADjjRgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUcQIAAAwihgBAABGESMAAMAoYgQAABhFjAAAAKOIEQAAYBQxAgAAjCJGAACAUQnFSEtLi0pKSuRyueT1erV3794xx27dulVz587Vd77zHeXn56uyslJvvPFGwhMGAACZxXaMtLW1qaGhQStXrlQgEFB1dbVqamrU3d0dc/yePXs0d+5ctbe3q6OjQzfffLNqa2sVCATOe/IAACD9OSzLsuzsMGfOHJWXl6u1tXV0W2lpqRYsWKDm5ua4vsf3vvc91dXV6YknnohrfDgcltvtVigUUn5+vp3pAgAAQ+L9/W3rysipU6fU0dEhn88Xsd3n82nfvn1xfY/h4WENDAxo8uTJY44ZHBxUOByOeAEAgMxkK0b6+/s1NDQkj8cTsd3j8ai3tzeu7/G73/1OJ06c0J133jnmmObmZrnd7tFXUVGRnWkCAIA0ktANrA6HI+Jry7KitsXy8ssv61e/+pXa2to0bdq0Mcc1NTUpFAqNvnp6ehKZJgAASAM5dgZPnTpV2dnZUVdB+vr6oq6WfFNbW5uWLFmiV155RT/5yU/OOtbpdMrpdNqZGgAASFO2rozk5ubK6/XK7/dHbPf7/aqqqhpzv5dffln33nuv/vKXv+j2229PbKYAACAj2boyIkmNjY1auHChKioqVFlZqQ0bNqi7u1v19fWSzrzF8tlnn+nFF1+UdCZEFi1apD/84Q+64YYbRq+qXHTRRXK73eO4FAAAkI5sx0hdXZ2OHTumNWvWKBgMqqysTO3t7SouLpYkBYPBiGeOPPvsszp9+rQefPBBPfjgg6PbFy9erBdeeOH8VwAAANKa7eeMmMBzRgAASD9Jec4IAADAeCNGAACAUcQIAAAwihgBAABG2f40DQBg4g0NWzrQdVx9Ayc1Lc+l2SWTlZ117idfA+mAGAGAFLfjcFCrt3cqGDo5uq3Q7dKq2lmaV1ZocGbA+OBtGgBIYTsOB7X0pXcjQkSSekMntfSld7XjcNDQzIDxQ4wAQIoaGra0enunYj0MamTb6u2dGhpO+cdFAWdFjABAijrQdTzqisjXWZKCoZM60HV84iYFJAExAgApqm9g7BBJZByQqogRAEhR0/Jc4zoOSFXECACkqNklk1XodmmsD/A6dOZTNbNLJk/ktIBxR4wAQIrKznJoVe0sSYoKkpGvV9XO4nkjSHvECACksHllhWq9p1wF7si3YgrcLrXeU85zRpAReOgZAKS4eWWFmjurgCewImMRIwCQBrKzHKq8aorpaQBJwds0AADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMIoYAQAARhEjAADAqBzTEwCAZBoatnSg67j6Bk5qWp5Ls0smKzvLYXpaAL6GGAGQsXYcDmr19k4FQydHtxW6XVpVO0vzygoNzgzA1/E2DYCMtONwUEtfejciRCSpN3RSS196VzsOBw3NDMA3ESMAMs7QsKXV2ztlxfi3kW2rt3dqaDjWCAATjRgBkHEOdB2PuiLydZakYOikDnQdn7hJARgTMQIg4/QNjB0iiYwDkFwX7A2smXSHfaashXWklnRex7Q817iOA5BcF2SMZNId9pmyFtaRWtJ9HbNLJqvQ7VJv6GTM+0YckgrcZwILgHkJvU3T0tKikpISuVwueb1e7d2796zjd+/eLa/XK5fLpSuvvFLr169PaLLjIZPusM+UtbCO1JIJ68jOcmhV7SxJZ8Lj60a+XlU7K22u9ACZznaMtLW1qaGhQStXrlQgEFB1dbVqamrU3d0dc3xXV5fmz5+v6upqBQIBPf7441q2bJm2bNly3pO3K5PusM+UtbCO1JIp65CkeWWFar2nXAXuyLdiCtwutd5TnhZXeIALhe23adauXaslS5bo/vvvlyStW7dOb7zxhlpbW9Xc3Bw1fv369br88su1bt06SVJpaakOHjyop556Snfcccf5zd4mO3fYV141ZeImloBMWQvrSC2Zso4R88oKNXdWQdre+wJcKGzFyKlTp9TR0aEVK1ZEbPf5fNq3b1/Mffbv3y+fzxex7bbbbtPGjRv11VdfadKkSVH7DA4OanBwcPTrcDhsZ5pjyqQ77DNlLawjtWTKOr4uO8uRFuEEXMhsvU3T39+voaEheTyeiO0ej0e9vb0x9+nt7Y05/vTp0+rv74+5T3Nzs9xu9+irqKjIzjTHlEl32GfKWlhHasmUdQBILwndwOpwRF7itCwratu5xsfaPqKpqUmhUGj01dPTk8g0o4zcYT/WTB0684mBdLjDPlPWwjpSS6asA0B6sRUjU6dOVXZ2dtRVkL6+vqirHyMKCgpijs/JydGUKbEvnTqdTuXn50e8xkMm3WGfKWthHaklU9YBIL3YipHc3Fx5vV75/f6I7X6/X1VVVTH3qaysjBq/c+dOVVRUxLxfJNky6Q77TFkL60gtmbIOAOnDYY28ZxKntrY2LVy4UOvXr1dlZaU2bNigP//5z3rvvfdUXFyspqYmffbZZ3rxxRclnflob1lZmR544AH97Gc/0/79+1VfX6+XX3457k/ThMNhud1uhUKhcbtKks5Pl/ymTFkL60gtmbIOAObE+/vb9kd76+rqdOzYMa1Zs0bBYFBlZWVqb29XcXGxJCkYDEY8c6SkpETt7e1avny5nnnmGU2fPl1//OMfJ/xjvd+USXfYZ8paWEdqyZR1AEh9tq+MmJCMKyMAACC54v39zV/tBQAARhEjAADAKGIEAAAYRYwAAACjiBEAAGAUMQIAAIwiRgAAgFHECAAAMMr2E1hNGHkuWzgcNjwTAAAQr5Hf2+d6vmpaxMjAwIAkqaioyPBMAACAXQMDA3K73WP+e1o8Dn54eFhHjx5VXl6eHI7x+0Nd4XBYRUVF6unp4THz44xjmzwc2+Th2CYHxzV5Uv3YWpalgYEBTZ8+XVlZY98ZkhZXRrKysjRjxoykff/8/PyU/I+YCTi2ycOxTR6ObXJwXJMnlY/t2a6IjOAGVgAAYBQxAgAAjLqgY8TpdGrVqlVyOp2mp5JxOLbJw7FNHo5tcnBckydTjm1a3MAKAAAy1wV9ZQQAAJhHjAAAAKOIEQAAYBQxAgAAjMr4GGlpaVFJSYlcLpe8Xq/27t171vG7d++W1+uVy+XSlVdeqfXr10/QTNOPnWO7a9cuORyOqNe///3vCZxx6tuzZ49qa2s1ffp0ORwOvfbaa+fch3M2PnaPLedsfJqbm/XDH/5QeXl5mjZtmhYsWKAPPvjgnPtx3p5bIsc2Xc/bjI6RtrY2NTQ0aOXKlQoEAqqurlZNTY26u7tjju/q6tL8+fNVXV2tQCCgxx9/XMuWLdOWLVsmeOapz+6xHfHBBx8oGAyOvq6++uoJmnF6OHHihK677jo9/fTTcY3nnI2f3WM7gnP27Hbv3q0HH3xQ77zzjvx+v06fPi2fz6cTJ06MuQ/nbXwSObYj0u68tTLY7Nmzrfr6+oht1157rbVixYqY4x977DHr2muvjdj2wAMPWDfccEPS5piu7B7bN99805Jk/fe//52A2WUGSda2bdvOOoZzNjHxHFvO2cT09fVZkqzdu3ePOYbzNjHxHNt0PW8z9srIqVOn1NHRIZ/PF7Hd5/Np3759MffZv39/1PjbbrtNBw8e1FdffZW0uaabRI7tiOuvv16FhYW69dZb9eabbyZzmhcEztnk45y1JxQKSZImT5485hjO28TEc2xHpNt5m7Ex0t/fr6GhIXk8nojtHo9Hvb29Mffp7e2NOf706dPq7+9P2lzTTSLHtrCwUBs2bNCWLVu0detWzZw5U7feeqv27NkzEVPOWJyzycM5a59lWWpsbNRNN92ksrKyMcdx3toX77FN1/M2Lf5q7/lwOBwRX1uWFbXtXONjbYe9Yztz5kzNnDlz9OvKykr19PToqaee0o9+9KOkzjPTcc4mB+esfQ899JD++c9/6q233jrnWM5be+I9tul63mbslZGpU6cqOzs76v/U+/r6oop8REFBQczxOTk5mjJlStLmmm4SObax3HDDDfroo4/Ge3oXFM7ZicU5O7aHH35Yr7/+ut58803NmDHjrGM5b+2xc2xjSYfzNmNjJDc3V16vV36/P2K73+9XVVVVzH0qKyujxu/cuVMVFRWaNGlS0uaabhI5trEEAgEVFhaO9/QuKJyzE4tzNpplWXrooYe0detW/f3vf1dJSck59+G8jU8ixzaWtDhvjd06OwH++te/WpMmTbI2btxodXZ2Wg0NDdYll1xiffzxx5ZlWdaKFSushQsXjo4/cuSIdfHFF1vLly+3Ojs7rY0bN1qTJk2yXn31VVNLSFl2j+3vf/97a9u2bdaHH35oHT582FqxYoUlydqyZYupJaSkgYEBKxAIWIFAwJJkrV271goEAtYnn3xiWRbn7Pmwe2w5Z+OzdOlSy+12W7t27bKCweDo68svvxwdw3mbmESObbqetxkdI5ZlWc8884xVXFxs5ebmWuXl5REfiVq8eLH14x//OGL8rl27rOuvv97Kzc21rrjiCqu1tXWCZ5w+7BzbX//619ZVV11luVwu69vf/rZ10003WX/7298MzDq1jXws75uvxYsXW5bFOXs+7B5bztn4xDqmkqznn39+dAznbWISObbpet46LOv/7xoCAAAwIGPvGQEAAOmBGAEAAEYRIwAAwChiBAAAGEWMAAAAo4gRAABgFDECAACMIkYAAIBRxAgAADCKGAEAAEYRIwAAwChiBAAAGPV/WwzCnc8KGdMAAAAASUVORK5CYII=",
      "text/plain": [
       "<Figure size 640x480 with 1 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.scatter(im, number_of_collapses/number_of_analyses)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 99,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "teta                beta                log_likelihood_sum  \n"
     ]
    }
   ],
   "source": [
    "#https://stackoverflow.com/questions/21610198/runtimewarning-divide-by-zero-encountered-in-log\n",
    "np.seterr(divide = 'ignore') #equal to IFERROR function of EXCEL\n",
    "im_log = np.log(im)\n",
    "\n",
    "FORMAT_STRING = \"{:<20}{:<20}{:<20}\"\n",
    "print(FORMAT_STRING.format(\"teta\", \"beta\", \"log_likelihood_sum\"))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Implementation to avoid underflow or overflow"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 100,
   "metadata": {},
   "outputs": [],
   "source": [
    "def neg_log_likelihood_sum(params, im_l, no_a, no_c):\n",
    "    \n",
    "    teta, beta = params\n",
    "\n",
    "    # log_theoretical_fragility_function = stats.norm.logcdf(im_l, np.log(teta), beta) #fisrt approach\n",
    "    # theoretical_fragility_function = np.exp(log_theoretical_fragility_function) #first approach\n",
    "    log_theoretical_fragility_function = stats.lognorm.logcdf(x=im, s=beta, loc=0, scale=teta) #second approach\n",
    "    theoretical_fragility_function = np.exp(log_theoretical_fragility_function) #second approach\n",
    "\n",
    "    log_likelihood = stats.binom.logpmf(no_c, no_a, theoretical_fragility_function)\n",
    "    log_likelihood_sum = np.sum(log_likelihood)\n",
    "\n",
    "    print(FORMAT_STRING.format(teta, beta, log_likelihood_sum))\n",
    "    \n",
    "    return -log_likelihood_sum"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 101,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "1.6                 0.4                 -185333.4207731888  \n",
      "1.6800000000000002  0.4                 -173166.32571673577 \n",
      "1.6                 0.42000000000000004 -186782.70388029766 \n",
      "1.6800000000000002  0.38                -171610.21971362675 \n",
      "1.7199999999999998  0.3600000000000001  -inf                \n",
      "1.7600000000000002  0.38                -164022.55083630758 \n",
      "1.8400000000000007  0.3699999999999999  -inf                \n",
      "1.7600000000000002  0.36                -inf                \n",
      "1.7000000000000002  0.39                -170010.7163315206  \n",
      "1.7800000000000002  0.39                -163780.74273908712 \n",
      "1.83                0.395               -162604.58593535333 \n",
      "1.8900000000000001  0.385               -161644.34413598548 \n",
      "1.9850000000000003  0.38250000000000006 -165597.8888591954  \n",
      "1.96                0.4                 -165358.269011201   \n",
      "1.81                0.385               -162101.84610934078 \n",
      "1.87                0.375               -inf                \n",
      "1.84                0.39                -161948.1200640003  \n",
      "1.9200000000000004  0.39                -162832.68457058232 \n",
      "1.8375000000000001  0.38625             -161622.1194245237  \n",
      "1.8875              0.38125             -161293.66918416545 \n",
      "1.9112500000000003  0.37687499999999985 -161546.29959350664 \n",
      "1.835               0.38249999999999995 -161317.9013317809  \n",
      "1.885               0.37749999999999995 -160963.8379871396  \n",
      "1.90875             0.37312499999999993 -inf                \n",
      "1.9375              0.37625             -inf                \n",
      "1.860625            0.38093749999999993 -161018.38276114786 \n",
      "1.858125            0.3771874999999999  -160714.72435304403 \n",
      "1.8434374999999998  0.37515624999999986 -160600.40609730277 \n",
      "1.8678125           0.3717187499999999  -inf                \n",
      "1.862421875         0.3786328124999999  -160834.8837568218  \n",
      "1.8208593749999997  0.3762890624999998  -161064.64930416452 \n",
      "1.86896484375       0.3771972656249999  -160754.68734608503 \n",
      "1.84998046875       0.37372070312499983 -inf                \n",
      "1.8593115234375     0.37740478515624987 -160732.4292632578  \n",
      "1.8337841796874998  0.37536376953124984 -160746.69632007775 \n",
      "1.842579345703125   0.3758221435546874  -160677.58206215495 \n",
      "1.8267053222656249  0.37357360839843734 -inf                \n",
      "1.8511599731445312  0.37644699096679674 -160669.83053328106 \n",
      "1.852018127441406   0.37578109741210913 -160616.2912829974  \n",
      "1.8442956542968747  0.3744903564453123  -inf                \n",
      "1.8494438934326172  0.37595783233642566 -160638.64386846585 \n",
      "1.8460117340087887  0.3749795150756834  -inf                \n",
      "1.84858585357666    0.3757132530212401  -160627.31563460286 \n",
      "1.846869773864746   0.37522409439086896 -160585.04005901638 \n",
      "1.8460117340087887  0.3749795150756833  -inf                \n",
      "1.8382891464233397  0.3745992469787596  -inf                \n",
      "1.8485858821868895  0.37548563480377173 -160611.43578731988 \n",
      "1.8417213916778563  0.37489470958709703 -inf                \n",
      "1.846869759559631   0.37533790349960305 -160596.64974095093 \n",
      "1.8503020334243772  0.3754057478904722  -160588.6345872858  \n",
      "1.850302047729492   0.37529193878173805 -160577.05804199737 \n",
      "1.8520181918144223  0.37526895642280556 -inf                \n",
      "1.8468697881698608  0.37511028528213486 -inf                \n",
      "1.849443972110748   0.37533188223838787 -160598.14607515908 \n",
      "1.848585910797119   0.37525801658630353 -160580.12917045073 \n",
      "1.8503020405769346  0.37534884333610513 -160596.8357951968  \n",
      "1.8485859179496764  0.37520111203193646 -inf                \n",
      "1.84987300992012    0.37531191051006296 -160587.54387773256 \n",
      "1.849014948606491   0.3752380448579786  -160569.29637406522 \n",
      "1.8485859179496762  0.3752011120319365  -inf                \n",
      "1.8507310855388641  0.37527196705341315 -160566.6885092692  \n",
      "1.8518036729097362  0.375278942286968   -160561.04637763344 \n",
      "1.8505165737867353  0.3752250483632086  -inf                \n",
      "1.850355679243803   0.3752752161771057  -160570.81079978147 \n",
      "1.8504629422724244  0.37524177096784095 -inf                \n",
      "1.8503824950009582  0.3752668548747895  -160567.6878454412  \n",
      "1.8531712193042034  0.3753077523037789  -160561.81010198244 \n",
      "1.8545923972129814  0.3753198397159574  -160557.6318241721  \n",
      "1.8566973483189928  0.37534633213654134 -160556.7223514242  \n",
      "1.8553298019245257  0.37531752211973046 -inf                \n",
      "1.853710864959284   0.3753101947577668  -160559.33020541543 \n",
      "1.8586045403685407  0.37537758460734016 -160560.77095004587 \n",
      "1.8569043235038394  0.3753529240272472  -160558.16395770665 \n",
      "1.8598908068635482  0.37538906140602174 -160561.43350676697 \n",
      "1.85525585043535    0.3753299114198305  -160557.64740385284 \n",
      "1.8550488752505037  0.3753233195291247  -160556.41569848734 \n",
      "1.8541211511238354  0.3753085172800634  -160556.33852269995 \n",
      "1.8555626490074781  0.37532493799677424 -160554.40356490976 \n",
      "1.8557160482935422  0.3753224512852462  -inf                \n",
      "1.8529864518123205  0.37528712314029633 -160555.82562593117 \n",
      "1.854427949695963   0.37530354385700715 -inf                \n",
      "1.8541978507668673  0.37530727392429936 -160555.46802122015 \n",
      "1.8567740479620252  0.3753450887807772  -160555.9744587567  \n",
      "1.8539333508497466  0.37530161455041655 -160555.0311427349  \n",
      "1.8552981490903577  0.37531927862289144 -160553.74365974223 \n",
      "1.855848298252103   0.37532528097218754 -inf                \n",
      "1.8569274472480894  0.3753426020692491  -160554.4896440976  \n",
      "1.8561789231485035  0.37533235518954094 -160554.10613721653 \n",
      "1.8559144232313833  0.37532669581565814 -160553.34570065408 \n",
      "1.8560903103433364  0.37532757472510014 -inf                \n",
      "1.8550336491732373  0.3753136192490087  -inf                \n",
      "1.855892604654687   0.3753276712044079  -160553.78546902555 \n",
      "1.8553199676670538  0.3753183032341417  -inf                \n",
      "1.8557494454077788  0.37532532921184136 -160553.64408914855 \n",
      "1.8563657195488044  0.37533274640460806 -inf                \n",
      "1.8555650417049694  0.3753226455683206  -160553.5962240324  \n",
      "1.855730019528574   0.3753240121721374  -inf                \n",
      "1.8557445889379776  0.3753249999519154  -160553.55280571483 \n",
      "1.8560939704643913  0.3753290501992529  -inf                \n",
      "1.8556972738948248  0.3753242467260537  -160553.51622433076 \n",
      "1.8558671081882303  0.37532594258979646 -inf                \n",
      "1.8557752187505407  0.37532523561138564 -160553.48977347527 \n",
      "1.8559923680870996  0.3753276847009901  -inf                \n",
      "1.8557710474428935  0.37532510621978776 -160553.46455150883 \n",
      "1.8559102519237358  0.37532656642406026 -inf                \n",
      "1.8558089770438395  0.3753255683145543  -160553.4463039911  \n",
      "1.8559523528323294  0.3753271579104247  -inf                \n",
      "1.8558163737902524  0.375325619142447   -160553.42921306822 \n",
      "1.855921819977796   0.3753267466435508  -160553.32909185917 \n",
      "1.8559782414447743  0.3753273358080491  -inf                \n",
      "1.856019869418927   0.37532782331676195 -inf                \n",
      "1.8558672476974212  0.37532617018602576 -160553.38250313845 \n",
      "1.8559689955117582  0.3753272722731832  -inf                \n",
      "1.8558926846510055  0.37532644570781515 -160553.35974773992 \n",
      "1.8559435585581738  0.3753269967513938  -inf                \n",
      "1.8559054031277977  0.37532658346870984 -160553.34851992398 \n",
      "1.8559308400813816  0.3753268589904991  -inf                \n",
      "1.8559117623661936  0.3753266523491572  -160553.34294348577 \n",
      "1.855919159112606   0.37532670317704986 -inf                \n",
      "1.855915607201689   0.37532669765600607 -160553.34085482228 \n",
      "1.8559256648132911  0.37532679195039975 -inf                \n",
      "1.855915237977968   0.3753266872494678  -160553.33895319005 \n",
      "1.8559214507540749  0.37532673623701257 -inf                \n",
      "1.8559170680897854  0.3753267073012577  -160553.33743634005 \n",
      "1.8559236500896135  0.3753267666953407  -inf                \n",
      "1.8559173410058794  0.3753267071109361  -160553.33610640588 \n",
      "1.8559220928938902  0.37532674645322917 -inf                \n",
      "1.8559183242908115  0.3753267170892506  -160553.3350164735  \n",
      "1.855922803262728   0.3753267566218653  -inf                \n",
      "1.8559187065700915  0.37532671948866836 -160553.3340792378  \n",
      "1.855922202257076   0.3753267490429686  -160553.32815544982 \n",
      "1.8559241412402079  0.37532676501982754 -inf                \n",
      "1.8559253156647804  0.37532677619785104 -inf                \n",
      "1.8559203588437638  0.37532673366596403 -160553.3313505968  \n",
      "1.855923663391108   0.37532676202055537 -inf                \n",
      "1.8559211849805999  0.37532674075461187 -160553.32998690908 \n",
      "1.855922837254272   0.37532675493190754 -inf                \n",
      "1.8559215980490178  0.37532674429893575 -160553.32930522345 \n",
      "1.855922424185854   0.37532675138758365 -inf                \n",
      "1.8559218045832269  0.3753267460710977  -160553.32896442013 \n",
      "Optimization terminated successfully.\n",
      "         Current function value: 160553.328155\n",
      "         Iterations: 71\n",
      "         Function evaluations: 140\n",
      " final_simplex: (array([[1.8559222 , 0.37532675],\n",
      "       [1.8559218 , 0.37532675],\n",
      "       [1.85592182, 0.37532675]]), array([160553.32815545, 160553.32896442, 160553.32909186]))\n",
      "           fun: 160553.32815544982\n",
      "       message: 'Optimization terminated successfully.'\n",
      "          nfev: 140\n",
      "           nit: 71\n",
      "        status: 0\n",
      "       success: True\n",
      "             x: array([1.8559222 , 0.37532675])\n"
     ]
    }
   ],
   "source": [
    "guess_values = (1.6, 0.4) #guess_values = (1.6, 0.1) for Baker, = (1.6, 0.4) for my data, while guess_values = (8, 0.1) for Qeshta et al. (2021) data. Be aware to use 'reasonable' guess data.\n",
    "params_bounds = ((0.01, None), (0.01, None)) #teta and beta must be positive because of the problem physics\n",
    "desired_tolerance = 10**-3\n",
    "number_of_iterations = 1000\n",
    "\n",
    "# https://www.udacity.com/blog/2020/12/how-to-create-partial-functions-in-python.html#:~:text=What%20is%20a%20Partial%20Function,extension%20of%20another%20specified%20function.\n",
    "neg_log_likelihood_sum_partial = partial(neg_log_likelihood_sum, im_l=im_log, no_a=number_of_analyses, no_c=number_of_collapses)\n",
    "\n",
    "#https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize.html\n",
    "res = optimize.minimize(neg_log_likelihood_sum_partial, x0 = guess_values, \n",
    "                method=\"Nelder-Mead\", options={'maxiter':number_of_iterations, \n",
    "                'disp':True}, bounds = params_bounds, tol = desired_tolerance)\n",
    "print(res)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Try to including 0"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 82,
   "metadata": {},
   "outputs": [],
   "source": [
    "xx=1**-11\n",
    "beta=0.4\n",
    "teta=1.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 83,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.15537184730329506\n"
     ]
    }
   ],
   "source": [
    "theoretical_fragility_function = stats.norm(np.log(teta), beta).cdf(np.log(xx))\n",
    "print(theoretical_fragility_function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 84,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.8619340202450088\n",
      "0.155371847303295\n"
     ]
    }
   ],
   "source": [
    "theoretical_fragility_function = stats.norm.logcdf(np.log(xx), np.log(teta), beta)\n",
    "print(theoretical_fragility_function)\n",
    "print(np.exp(theoretical_fragility_function))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 85,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.155371847303295\n"
     ]
    }
   ],
   "source": [
    "#https://docs.scipy.org/doc/scipy/reference/generated/scipy.stats.lognorm.html\n",
    "theoretical_fragility_function = stats.lognorm.cdf(x=xx, s=beta, loc=0, scale=teta)\n",
    "print(theoretical_fragility_function)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 86,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-1.861934020245009\n",
      "0.15537184730329498\n"
     ]
    }
   ],
   "source": [
    "theoretical_fragility_function = stats.lognorm.logcdf(x=xx, s=beta, loc=0, scale=teta)\n",
    "print(theoretical_fragility_function)\n",
    "print(np.exp(theoretical_fragility_function))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 87,
   "metadata": {},
   "outputs": [],
   "source": [
    "no_c = 69\n",
    "no_a = 100_000"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 88,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-16434.52479231295\n",
      "0.0\n"
     ]
    }
   ],
   "source": [
    "log_likelihood = stats.binom.logpmf(no_c, no_a, np.exp(theoretical_fragility_function))\n",
    "print(log_likelihood)\n",
    "print(np.exp(log_likelihood))"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "base",
   "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.13"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "ad2bdc8ecc057115af97d19610ffacc2b4e99fae6737bb82f5d7fb13d2f2c186"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}