wasting_transitions/alibow_vicious_cycle/symbolic_prevalence_equation_solver.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "024b7d89",
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np, pandas as pd\n",
    "import sympy as sym\n",
    "from sympy import symbols, Matrix, solve, simplify\n",
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "f50255c0",
   "metadata": {},
   "outputs": [],
   "source": [
    "#### initialize unknown variable symbols (30 in total)\n",
    "\n",
    "# susceptible to diarrhea prevalence variables\n",
    "p_S1, p_S2, p_S3, p_S4 = symbols('p_S1 p_S2 p_S3 p_S4')\n",
    "\n",
    "# infected with diarrhea prevalence variables\n",
    "p_D1, p_D2, p_D3, p_D4 = symbols('p_D1 p_D2 p_D3 p_D4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "718f90d2",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "unknowns = [p_S1, p_S2, p_S3, \n",
    "           p_D1, p_D2, p_D3, p_S4, p_D4\n",
    "]\n",
    "len(unknowns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "0faddfb3",
   "metadata": {},
   "outputs": [],
   "source": [
    "### initialize known variable symbols\n",
    "\n",
    "# diarrhea cause model variables\n",
    "prevalence_c302 = symbols('prevalence_c302') \n",
    "\n",
    "# wasting risk exposure values\n",
    "p_4, p_3, p_2, p_1 = symbols('p_4 p_3 p_2 p_1')\n",
    "\n",
    "# wasting exposure prevalence ratios by diarrheal state\n",
    "PR_1, PR_2, PR_3, PR_4 = symbols('PR_1 PR_2 PR_3 PR_4')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "344b3a40",
   "metadata": {},
   "outputs": [],
   "source": [
    "#[{unknown:coefficient, unknown:coefficient}, right hand side]\n",
    "\n",
    "eq1 = [{p_D1:1/prevalence_c302, p_S1:-PR_1/(1-prevalence_c302)}, 0]\n",
    "eq2 = [{p_D2:1/prevalence_c302, p_S2:-PR_2/(1-prevalence_c302)}, 0]\n",
    "eq3 = [{p_D3:1/prevalence_c302, p_S3:-PR_3/(1-prevalence_c302)}, 0]\n",
    "eq4 = [{p_D4:1/prevalence_c302, p_S4:-PR_4/(1-prevalence_c302)}, 0]\n",
    "eq5 = [{p_D1:1, p_D2:1, p_D3:1, p_D4:1}, prevalence_c302]\n",
    "eq6 = [{p_S1:1, p_D1:1}, p_1]\n",
    "eq7 = [{p_S2:1, p_D2:1}, p_2]\n",
    "eq8 = [{p_S3:1, p_D3:1}, p_3]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1fc3e502",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "8"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eqs = [eq1,eq2,eq3,eq4,eq5,eq6,eq7,eq8]\n",
    "len(eqs)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "311288cf",
   "metadata": {},
   "outputs": [],
   "source": [
    "def build_matrix(eqns, unknowns):\n",
    "    \"\"\"\n",
    "    INPUT\n",
    "    ----\n",
    "    eqns: a list of equations. These are of the form:\n",
    "          [{x_i:a, x_j:b}, y] for the equation:\n",
    "          ax_i + bx_j = y\n",
    "    unknowns: a list of sympy unknowns\n",
    "    ----\n",
    "    OUTPUT\n",
    "    ----\n",
    "    A:  a matrix containing the coefficients of left-hand side of all eq in eqns.\n",
    "        - nrows = number of equations\n",
    "        - rcols = number of unknowns\n",
    "    b: an nx1 matrix containing the right-hand side of all the eqns\n",
    "    x: a sympy matrix of the unknowns\n",
    "    \n",
    "    Note that Ax = b\n",
    "    \"\"\"\n",
    "    n_eqns = len(eqns)\n",
    "    n_unknowns = len(unknowns)\n",
    "\n",
    "    # frame for matrix/LHS equation coefficients.\n",
    "    # nrows = n_eqns, ncols = n_unknowns\n",
    "    A = pd.DataFrame(\n",
    "        index = range(n_eqns),\n",
    "        columns = unknowns,\n",
    "        data = np.zeros([n_eqns,n_unknowns])\n",
    "    )\n",
    "\n",
    "    # frame for RHS of equations\n",
    "    b = pd.DataFrame(index = range(n_eqns), columns = ['val'])\n",
    "    \n",
    "    # define helper fn\n",
    "    def add_eq(terms, y, i):\n",
    "        \"\"\"\n",
    "        INPUT\n",
    "        ----\n",
    "        To input i^(th) equation y = ax_i + bx_j, add_eq wants:\n",
    "            - terms = {x_i:a, x_j:b}\n",
    "            - y = y\n",
    "            - i\n",
    "        ----\n",
    "        FUNCTION\n",
    "        ----\n",
    "        - Adds coefficients to matrix A.\n",
    "        - Adds y to row i of vector b\n",
    "        ----\n",
    "        EXAMPLE\n",
    "        ----\n",
    "        For y = ax_i + bx_j:\n",
    "        - adds \"a\" to column \"x_i\" row i of A\n",
    "        - adds \"b\" to column \"x_j\" row i of A\n",
    "        - adds \"y\" to row i of b\n",
    "        \"\"\"\n",
    "        for x in terms.keys():\n",
    "            A[x][i] = terms[x]\n",
    "        b.iloc[i] = y\n",
    "\n",
    "    # populate LHS/RHS\n",
    "    i = 0\n",
    "    for eq in eqns:\n",
    "        add_eq(terms=eq[0], y=eq[1], i=i)\n",
    "        i += 1\n",
    "\n",
    "    # convert to sympy matrices\n",
    "    A = sym.Matrix(A)\n",
    "    b = sym.Matrix(b)\n",
    "    x = sym.Matrix(unknowns) #vars to solve for\n",
    "\n",
    "    return A, x, b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "dde95361",
   "metadata": {},
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "/ihme/code/alibow/miniconda/envs/vivarium_ciff_sam/lib/python3.6/site-packages/ipykernel_launcher.py:57: SettingWithCopyWarning: \n",
      "A value is trying to be set on a copy of a slice from a DataFrame\n",
      "\n",
      "See the caveats in the documentation: https://pandas.pydata.org/pandas-docs/stable/user_guide/indexing.html#returning-a-view-versus-a-copy\n"
     ]
    }
   ],
   "source": [
    "A1, x1, b1 = build_matrix(eqs, unknowns)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "73a88984",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}- \\frac{PR_{1}}{1 - prevalence_{c302}} & 0.0 & 0.0 & \\frac{1}{prevalence_{c302}} & 0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & - \\frac{PR_{2}}{1 - prevalence_{c302}} & 0.0 & 0.0 & \\frac{1}{prevalence_{c302}} & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & - \\frac{PR_{3}}{1 - prevalence_{c302}} & 0.0 & 0.0 & \\frac{1}{prevalence_{c302}} & 0.0 & 0.0\\\\0.0 & 0.0 & 0.0 & 0.0 & 0.0 & 0.0 & - \\frac{PR_{4}}{1 - prevalence_{c302}} & \\frac{1}{prevalence_{c302}}\\\\0.0 & 0.0 & 0.0 & 1 & 1 & 1 & 0.0 & 1\\\\1 & 0.0 & 0.0 & 1 & 0.0 & 0.0 & 0.0 & 0.0\\\\0.0 & 1 & 0.0 & 0.0 & 1 & 0.0 & 0.0 & 0.0\\\\0.0 & 0.0 & 1 & 0.0 & 0.0 & 1 & 0.0 & 0.0\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[-PR_1/(1 - prevalence_c302),                         0.0,                         0.0, 1/prevalence_c302,               0.0,               0.0,                         0.0,               0.0],\n",
       "[                        0.0, -PR_2/(1 - prevalence_c302),                         0.0,               0.0, 1/prevalence_c302,               0.0,                         0.0,               0.0],\n",
       "[                        0.0,                         0.0, -PR_3/(1 - prevalence_c302),               0.0,               0.0, 1/prevalence_c302,                         0.0,               0.0],\n",
       "[                        0.0,                         0.0,                         0.0,               0.0,               0.0,               0.0, -PR_4/(1 - prevalence_c302), 1/prevalence_c302],\n",
       "[                        0.0,                         0.0,                         0.0,                 1,                 1,                 1,                         0.0,                 1],\n",
       "[                          1,                         0.0,                         0.0,                 1,               0.0,               0.0,                         0.0,               0.0],\n",
       "[                        0.0,                           1,                         0.0,               0.0,                 1,               0.0,                         0.0,               0.0],\n",
       "[                        0.0,                         0.0,                           1,               0.0,               0.0,                 1,                         0.0,               0.0]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "2922658a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 8)"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "A1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "22155584",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}p_{S1}\\\\p_{S2}\\\\p_{S3}\\\\p_{D1}\\\\p_{D2}\\\\p_{D3}\\\\p_{S4}\\\\p_{D4}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[p_S1],\n",
       "[p_S2],\n",
       "[p_S3],\n",
       "[p_D1],\n",
       "[p_D2],\n",
       "[p_D3],\n",
       "[p_S4],\n",
       "[p_D4]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "1eb5f591",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 1)"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "2404464d",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/latex": [
       "$\\displaystyle \\left[\\begin{matrix}0\\\\0\\\\0\\\\0\\\\prevalence_{c302}\\\\p_{1}\\\\p_{2}\\\\p_{3}\\end{matrix}\\right]$"
      ],
      "text/plain": [
       "Matrix([\n",
       "[              0],\n",
       "[              0],\n",
       "[              0],\n",
       "[              0],\n",
       "[prevalence_c302],\n",
       "[            p_1],\n",
       "[            p_2],\n",
       "[            p_3]])"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b1"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "51464709",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "(8, 1)"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "b1.shape"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "9630be7a",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "{p_S1: (-p_1*prevalence_c302 + p_1)/(PR_1*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_S2: (-p_2*prevalence_c302 + p_2)/(PR_2*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_S3: (-p_3*prevalence_c302 + p_3)/(PR_3*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_D1: PR_1*p_1*prevalence_c302/(PR_1*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_D2: PR_2*p_2*prevalence_c302/(PR_2*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_D3: PR_3*p_3*prevalence_c302/(PR_3*prevalence_c302 - prevalence_c302 + 1),\n",
       " p_S4: (PR_1*PR_2*PR_3*p_1*prevalence_c302**3 - PR_1*PR_2*PR_3*p_1*prevalence_c302**2 + PR_1*PR_2*PR_3*p_2*prevalence_c302**3 - PR_1*PR_2*PR_3*p_2*prevalence_c302**2 + PR_1*PR_2*PR_3*p_3*prevalence_c302**3 - PR_1*PR_2*PR_3*p_3*prevalence_c302**2 - PR_1*PR_2*PR_3*prevalence_c302**4 + PR_1*PR_2*PR_3*prevalence_c302**3 - PR_1*PR_2*p_1*prevalence_c302**3 + 2*PR_1*PR_2*p_1*prevalence_c302**2 - PR_1*PR_2*p_1*prevalence_c302 - PR_1*PR_2*p_2*prevalence_c302**3 + 2*PR_1*PR_2*p_2*prevalence_c302**2 - PR_1*PR_2*p_2*prevalence_c302 + PR_1*PR_2*prevalence_c302**4 - 2*PR_1*PR_2*prevalence_c302**3 + PR_1*PR_2*prevalence_c302**2 - PR_1*PR_3*p_1*prevalence_c302**3 + 2*PR_1*PR_3*p_1*prevalence_c302**2 - PR_1*PR_3*p_1*prevalence_c302 - PR_1*PR_3*p_3*prevalence_c302**3 + 2*PR_1*PR_3*p_3*prevalence_c302**2 - PR_1*PR_3*p_3*prevalence_c302 + PR_1*PR_3*prevalence_c302**4 - 2*PR_1*PR_3*prevalence_c302**3 + PR_1*PR_3*prevalence_c302**2 + PR_1*p_1*prevalence_c302**3 - 3*PR_1*p_1*prevalence_c302**2 + 3*PR_1*p_1*prevalence_c302 - PR_1*p_1 - PR_1*prevalence_c302**4 + 3*PR_1*prevalence_c302**3 - 3*PR_1*prevalence_c302**2 + PR_1*prevalence_c302 - PR_2*PR_3*p_2*prevalence_c302**3 + 2*PR_2*PR_3*p_2*prevalence_c302**2 - PR_2*PR_3*p_2*prevalence_c302 - PR_2*PR_3*p_3*prevalence_c302**3 + 2*PR_2*PR_3*p_3*prevalence_c302**2 - PR_2*PR_3*p_3*prevalence_c302 + PR_2*PR_3*prevalence_c302**4 - 2*PR_2*PR_3*prevalence_c302**3 + PR_2*PR_3*prevalence_c302**2 + PR_2*p_2*prevalence_c302**3 - 3*PR_2*p_2*prevalence_c302**2 + 3*PR_2*p_2*prevalence_c302 - PR_2*p_2 - PR_2*prevalence_c302**4 + 3*PR_2*prevalence_c302**3 - 3*PR_2*prevalence_c302**2 + PR_2*prevalence_c302 + PR_3*p_3*prevalence_c302**3 - 3*PR_3*p_3*prevalence_c302**2 + 3*PR_3*p_3*prevalence_c302 - PR_3*p_3 - PR_3*prevalence_c302**4 + 3*PR_3*prevalence_c302**3 - 3*PR_3*prevalence_c302**2 + PR_3*prevalence_c302 + prevalence_c302**4 - 4*prevalence_c302**3 + 6*prevalence_c302**2 - 4*prevalence_c302 + 1)/(PR_1*PR_2*PR_3*PR_4*prevalence_c302**3 - PR_1*PR_2*PR_4*prevalence_c302**3 + PR_1*PR_2*PR_4*prevalence_c302**2 - PR_1*PR_3*PR_4*prevalence_c302**3 + PR_1*PR_3*PR_4*prevalence_c302**2 + PR_1*PR_4*prevalence_c302**3 - 2*PR_1*PR_4*prevalence_c302**2 + PR_1*PR_4*prevalence_c302 - PR_2*PR_3*PR_4*prevalence_c302**3 + PR_2*PR_3*PR_4*prevalence_c302**2 + PR_2*PR_4*prevalence_c302**3 - 2*PR_2*PR_4*prevalence_c302**2 + PR_2*PR_4*prevalence_c302 + PR_3*PR_4*prevalence_c302**3 - 2*PR_3*PR_4*prevalence_c302**2 + PR_3*PR_4*prevalence_c302 - PR_4*prevalence_c302**3 + 3*PR_4*prevalence_c302**2 - 3*PR_4*prevalence_c302 + PR_4),\n",
       " p_D4: (-PR_1*PR_2*PR_3*p_1*prevalence_c302**3 - PR_1*PR_2*PR_3*p_2*prevalence_c302**3 - PR_1*PR_2*PR_3*p_3*prevalence_c302**3 + PR_1*PR_2*PR_3*prevalence_c302**4 + PR_1*PR_2*p_1*prevalence_c302**3 - PR_1*PR_2*p_1*prevalence_c302**2 + PR_1*PR_2*p_2*prevalence_c302**3 - PR_1*PR_2*p_2*prevalence_c302**2 - PR_1*PR_2*prevalence_c302**4 + PR_1*PR_2*prevalence_c302**3 + PR_1*PR_3*p_1*prevalence_c302**3 - PR_1*PR_3*p_1*prevalence_c302**2 + PR_1*PR_3*p_3*prevalence_c302**3 - PR_1*PR_3*p_3*prevalence_c302**2 - PR_1*PR_3*prevalence_c302**4 + PR_1*PR_3*prevalence_c302**3 - PR_1*p_1*prevalence_c302**3 + 2*PR_1*p_1*prevalence_c302**2 - PR_1*p_1*prevalence_c302 + PR_1*prevalence_c302**4 - 2*PR_1*prevalence_c302**3 + PR_1*prevalence_c302**2 + PR_2*PR_3*p_2*prevalence_c302**3 - PR_2*PR_3*p_2*prevalence_c302**2 + PR_2*PR_3*p_3*prevalence_c302**3 - PR_2*PR_3*p_3*prevalence_c302**2 - PR_2*PR_3*prevalence_c302**4 + PR_2*PR_3*prevalence_c302**3 - PR_2*p_2*prevalence_c302**3 + 2*PR_2*p_2*prevalence_c302**2 - PR_2*p_2*prevalence_c302 + PR_2*prevalence_c302**4 - 2*PR_2*prevalence_c302**3 + PR_2*prevalence_c302**2 - PR_3*p_3*prevalence_c302**3 + 2*PR_3*p_3*prevalence_c302**2 - PR_3*p_3*prevalence_c302 + PR_3*prevalence_c302**4 - 2*PR_3*prevalence_c302**3 + PR_3*prevalence_c302**2 - prevalence_c302**4 + 3*prevalence_c302**3 - 3*prevalence_c302**2 + prevalence_c302)/(PR_1*PR_2*PR_3*prevalence_c302**3 - PR_1*PR_2*prevalence_c302**3 + PR_1*PR_2*prevalence_c302**2 - PR_1*PR_3*prevalence_c302**3 + PR_1*PR_3*prevalence_c302**2 + PR_1*prevalence_c302**3 - 2*PR_1*prevalence_c302**2 + PR_1*prevalence_c302 - PR_2*PR_3*prevalence_c302**3 + PR_2*PR_3*prevalence_c302**2 + PR_2*prevalence_c302**3 - 2*PR_2*prevalence_c302**2 + PR_2*prevalence_c302 + PR_3*prevalence_c302**3 - 2*PR_3*prevalence_c302**2 + PR_3*prevalence_c302 - prevalence_c302**3 + 3*prevalence_c302**2 - 3*prevalence_c302 + 1)}"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "result = sym.solve(A1 * x1 - b1, x1)#, dict=True, set=True)\n",
    "result\n",
    "# NOTE: the p_S4 and p_D4 equations are complicated if they do not use\n",
    "# any of the other unknown variables, but can be simplified by using \n",
    "# those values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "7ac5dea0",
   "metadata": {},
   "outputs": [],
   "source": [
    "PR_1=1.041978\n",
    "PR_2=1.005639\n",
    "PR_3=1.006994\n",
    "PR_4=0.998314\n",
    "\n",
    "prevalence_c302 = 0.03\n",
    "p_1 = 0.016\n",
    "p_2 = 0.17\n",
    "p_3 = 0.4\n",
    "p_4 = 1-p_1-p_2-p_3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "a39c7569",
   "metadata": {},
   "outputs": [],
   "source": [
    "p_D1 = (PR_1 * p_1 * prevalence_c302) / (PR_1 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_D2 = (PR_2 * p_2 * prevalence_c302) / (PR_2 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_D3 = (PR_3 * p_3 * prevalence_c302) / (PR_3 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_S1 = (-p_1 * prevalence_c302 + p_1) / (PR_1 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_S2 = (-p_2 * prevalence_c302 + p_2) / (PR_2 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_S3 = (-p_3 * prevalence_c302 + p_3) / (PR_3 * prevalence_c302 - prevalence_c302 + 1)\n",
    "p_D4 = prevalence_c302 - p_D1 - p_D2 - p_D3\n",
    "p_S4 = (1 - prevalence_c302) - p_S1 - p_S2 - p_S3"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "b58c11c7",
   "metadata": {},
   "outputs": [],
   "source": [
    "import math"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "f7aff966",
   "metadata": {},
   "outputs": [],
   "source": [
    "assert math.isclose(p_D1+p_D2+p_D3+p_D4+p_S1+p_S2+p_S3+p_S4, 1, abs_tol=1e-8), 'Do not all sum to one'\n",
    "for i in [p_D1,p_D2,p_D3,p_D4,p_S1,p_S2,p_S3,p_S4]:\n",
    "    assert i > 0, f'{i} is negative'\n",
    "    assert i < 1, f'{i} is greater than one'"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "d9f574d5",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.0004995203740121915 0.00512789141460939 0.012081393082103513 0.012291195129274903\n"
     ]
    }
   ],
   "source": [
    "print(p_D1,\n",
    "     p_D2,\n",
    "     p_D3,\n",
    "     p_D4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "61de6f33",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "0.015500479625987808 0.16487210858539064 0.3879186069178965 0.401708804870725\n"
     ]
    }
   ],
   "source": [
    "print(p_S1,p_S2,p_S3,p_S4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "871a4f00",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.40970650430916344"
      ]
     },
     "execution_count": 24,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_D4 / (p_D1 + p_D2 + p_D3 + p_D4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "6c50a584",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.41413278852652063"
      ]
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "p_S4 / (p_S1 + p_S2 + p_S3 + p_S4)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "4da45555",
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "f96db554",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "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.6.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}