notebooks/conditional-probability-tables-with-constraints.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Generating conditional probability tables subject to constraints\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import os\n",
    "from pathlib import Path\n",
    "\n",
    "from itertools import product\n",
    "\n",
    "import numpy as np\n",
    "import pandas as pd\n",
    "\n",
    "from fake_data_for_learning.fake_data_for_learning import (\n",
    "    BayesianNodeRV, FakeDataBayesianNetwork, SampleValue\n",
    ")\n",
    "from fake_data_for_learning.utils import RandomCpt\n",
    "from fake_data_for_learning.probability_polytopes import (\n",
    "    MapMultidimIndexToLinear, ProbabilityPolytope, ExpectationConstraint\n",
    ")"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose we want to generate data from a discrete Bayesian network, such as\n",
    "\n",
    "Product -> Days <- Rating, \n",
    "\n",
    "where e.g. Product is the (insurance) product name, Rating is rating strength (i.e. market price / technical price) for a submission, and Days is the number of days to generate a quote for the submission.\n",
    "\n",
    "The number of entries in probability and conditional probability tables to define this Bayesian network is\n",
    "\n",
    "$ | Product | + | Rating | + | Product | \\times | Rating | \\times | Days |$.\n",
    "\n",
    "For example, let us define Industry and Rating as follows"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "product_values = ['financial', 'liability', 'property']\n",
    "product_type = BayesianNodeRV('product_type', np.array([0.2, 0.5, 0.3]), values=product_values)\n",
    "rating_values = range(2)\n",
    "rating = BayesianNodeRV('rating', np.array([0.3, 0.7]))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Suppose that Days is also discrete, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "days_values = range(4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Then if we choose the ordering of the conditional probability table axes as Product, Rating, Days, we can generate the entries of the conditional probability table for Days conditioned on Industry and Rating with `utils.RandomCpt`:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "random_cpt = RandomCpt(len(product_values), len(rating_values), len(days_values))\n",
    "X = random_cpt()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "0.9999999999999999"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "X[0, 0, :].sum()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "So the total number of probability table entries to specify is, as in the formula above,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "'Number of probability table entries: 29'"
      ]
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "f'Number of probability table entries: {len(product_values) + len(rating_values) +  (len(product_values) * len(rating_values) * len(days_values))}'"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "It would be nice to specify certain properties of the matrix without having to change entries individually. For example, we may want to insist that\n",
    "\n",
    "\\begin{equation*}\n",
    "E(D | P = property) = 3.5 \\\\\n",
    "E(D | P = financial) = 1.0 \\\\\n",
    "E(D | P= liability) = 2.0\n",
    "\\end{equation*}\n",
    "\n",
    "Denote the entries of the conditional probability table as \n",
    "\n",
    "$$(\\rho_{p, r | d})$$\n",
    "\n",
    "The the above constraints become\n",
    "\n",
    "\\begin{equation*}\n",
    "\\frac{1}{|R|} \\sum_{r, d} d \\, \\rho_{\\mathrm{property},\\, r\\, | d} = 3.5 \\\\\n",
    "\\frac{1}{|R|} \\sum_{r, d} d \\, \\rho_{\\mathrm{financial},\\, r\\, | d} = 1.0\\\\\n",
    "\\frac{1}{|R|} \\sum_{r, d} d \\, \\rho_{\\mathrm{liability},\\, r\\, | d} = 2.0.\n",
    "\\end{equation*}\n",
    "\n",
    "As $(\\rho)$ is a conditional probability table, we also have the constraints \n",
    "\n",
    "\\begin{equation*}\n",
    "0 \\leq \\rho_{p,\\,r\\,|d} \\leq 1 \\textrm{ for all }(p,\\,r,\\,d),\\\\\n",
    "\\sum_{d} \\rho_{p,\\,r,\\,| d} = 1 \\textrm{ for each pair } (p, \\, r)\n",
    "\\end{equation*}\n",
    "\n",
    "Together, these constraints define convex polytope contained in (probability) simplex $\\Delta_{R-1} \\subseteq \\mathbb{R}^{R}$, where $R = |Product | \\times | Rating | \\times | Days|$ (see e.g. Chapter 1 of *Lectures on Algebraic Statistics*, Drton, Sturmfels, Sullivant). This polytope is defined as an intersection of half-spaces, i.e. using the so-called *H-representation* of the polytope, see *Lectures on Polytopes* by Ziegler, Chapters 0 and 1.\n",
    "\n",
    "To generate a random (conditional) probability table to these constraints, the vertex-, or *V-representation* of the probability polytope $P$ is much more useful, because given the a vertex matrix $V$, where each column is a vertex of $P$ in $\\mathbb{R}^R$, and all points in $P$ can be obtained as\n",
    "\n",
    "$$\n",
    "\\begin{equation*}\n",
    "x = V \\cdot t\n",
    "\\end{equation*}\n",
    "$$\n",
    "\n",
    "where $t \\in \\mathbb{R}^N$, with $N$ being the number of vertices for $P$, and $t$ satisfying $0 \\leq t_i \\leq 1$, $\\sum t_i = 1$.\n",
    "\n",
    "Once we have determined the V-representation $V$, then the problem of generating conditional probability tables subject to our given expectation value constraints reduces to the much simpler problem of generating points on the non-negative quadrant of the unit (hyper) cube in $R^N$.\n",
    "\n",
    "Before we get to our goal of generating these probability tables for our hit ratio problem, let's look at elementary examples."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## (Conditional) Probability Polytopes\n",
    "\n",
    "The simplest example of a probability polytope is that of a Bernoulli random variable."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  1.]\n",
      " [-1. -1.]\n",
      " [ 1.  0.]\n",
      " [ 0.  1.]\n",
      " [-1. -0.]\n",
      " [-0. -1.]] \n",
      " [ 1. -1.  1.  1.  0.  0.]\n"
     ]
    }
   ],
   "source": [
    "bernoulli = ProbabilityPolytope(('outcome',), dict(outcome=range(2)))\n",
    "A, b = bernoulli.get_probability_half_planes()\n",
    "print(A, '\\n', b)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We convert the formulation A x <= b to the V-description"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0.],\n",
       "       [0., 1.]])"
      ]
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "bernoulli.get_vertex_representation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 0., 0.],\n",
       "       [0., 1., 0.],\n",
       "       [0., 0., 1.]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "tertiary = ProbabilityPolytope(('outcome',), dict(outcome=range(3)))\n",
    "tertiary.get_vertex_representation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1., 1., 0., 0.],\n",
       "       [0., 0., 1., 1.],\n",
       "       [1., 0., 0., 1.],\n",
       "       [0., 1., 1., 0.]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis = ProbabilityPolytope(\n",
    "    ('input', 'output'), dict(input=range(2), output=range(2))\n",
    ")\n",
    "conditional_bernoullis.get_vertex_representation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The benefit of having the vertex-representation (V-representation) of the probability polytope is that generating random (conditional) probability tables is straightforward, namely, we can get all elements of the probability polytope by taking combinations of the vertex (column) vectors.\n",
    "\n",
    "In the flattened coordinates, we have, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([0.27086489, 0.72913511, 0.74258253, 0.25741747])"
      ]
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis.generate_flat_random_cpt()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the multidimensional coordinates for conditional probability tables here, we have e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.58372782, 0.41627218],\n",
       "       [0.71527302, 0.28472698]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis.generate_random_cpt()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Adding contraints on conditional expectation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "conditional_bernoullis.set_expectation_constraints(\n",
    "    [ExpectationConstraint(equation=dict(input=1), moment=1, value=0.5)]\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[2, 3]"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis.get_expect_equations_col_indices(conditional_bernoullis.expect_constraints[0].equation)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[1. , 0. ],\n",
       "       [0. , 1. ],\n",
       "       [0.5, 0.5],\n",
       "       [0.5, 0.5]])"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis.get_vertex_representation()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.76396714, 0.23603286],\n",
       "       [0.5       , 0.5       ]])"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "conditional_bernoullis.generate_random_cpt()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[0.5, 0.5, 0.5, 0.5, 1. , 1. , 1. , 1. ],\n",
       "       [0.5, 0.5, 0.5, 0.5, 0. , 0. , 0. , 0. ],\n",
       "       [0. , 0. , 1. , 1. , 1. , 1. , 0. , 0. ],\n",
       "       [1. , 1. , 0. , 0. , 0. , 0. , 1. , 1. ],\n",
       "       [1. , 1. , 1. , 1. , 0.5, 0.5, 0.5, 0.5],\n",
       "       [0. , 0. , 0. , 0. , 0.5, 0.5, 0.5, 0.5],\n",
       "       [1. , 0. , 1. , 0. , 1. , 0. , 0. , 1. ],\n",
       "       [0. , 1. , 0. , 1. , 0. , 1. , 1. , 0. ]])"
      ]
     },
     "execution_count": 17,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "two_input_constrained_polytope = ProbabilityPolytope(\n",
    "    ('input', 'more_input', 'output'),\n",
    "    dict(input=['hi', 'low'], more_input=range(2), output=range(2))\n",
    ")\n",
    "two_input_constrained_polytope.set_expectation_constraints(\n",
    "    [ExpectationConstraint(equation=dict(more_input=0), moment=1, value=0.25)]\n",
    ")\n",
    "two_input_constrained_polytope.get_vertex_representation()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Hit rate polytope again"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[[0.88064391, 0.06169566, 0.03086844, 0.02679199],\n",
       "        [0.87308745, 0.07564295, 0.03326026, 0.01800934]],\n",
       "\n",
       "       [[0.40923928, 0.35176295, 0.14170279, 0.09729497],\n",
       "        [0.44115411, 0.32982671, 0.14393771, 0.08508147]],\n",
       "\n",
       "       [[0.66932198, 0.18715179, 0.08308226, 0.06044397],\n",
       "        [0.71299054, 0.15674574, 0.08218511, 0.0480786 ]]])"
      ]
     },
     "execution_count": 18,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "days_polytope = ProbabilityPolytope(\n",
    "    ('product', 'rating', 'days'),\n",
    "    coords = {\n",
    "        'product': product_values, \n",
    "        'rating': rating_values, \n",
    "        'days': days_values\n",
    "    }\n",
    ")\n",
    "days_polytope.set_expectation_constraints(\n",
    "    [\n",
    "        ExpectationConstraint(equation=dict(product='financial'), moment=1, value=0.2),\n",
    "        ExpectationConstraint(equation=dict(product='liability'), moment=1, value=0.9),\n",
    "        ExpectationConstraint(equation=dict(product='property'), moment=1, value=0.5),\n",
    "    ]\n",
    ")\n",
    "days_cpt = days_polytope.generate_random_cpt()\n",
    "days_cpt"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Now we create our Bayesian network with desired constraints on some expectation values"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<style scoped>\n",
       "    .dataframe tbody tr th:only-of-type {\n",
       "        vertical-align: middle;\n",
       "    }\n",
       "\n",
       "    .dataframe tbody tr th {\n",
       "        vertical-align: top;\n",
       "    }\n",
       "\n",
       "    .dataframe thead th {\n",
       "        text-align: right;\n",
       "    }\n",
       "</style>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th></th>\n",
       "      <th>product_type</th>\n",
       "      <th>rating</th>\n",
       "      <th>days</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>0</th>\n",
       "      <td>financial</td>\n",
       "      <td>1</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>1</th>\n",
       "      <td>property</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>2</th>\n",
       "      <td>financial</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>3</th>\n",
       "      <td>liability</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>4</th>\n",
       "      <td>property</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>5</th>\n",
       "      <td>property</td>\n",
       "      <td>0</td>\n",
       "      <td>1</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>6</th>\n",
       "      <td>liability</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>7</th>\n",
       "      <td>property</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>8</th>\n",
       "      <td>financial</td>\n",
       "      <td>0</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>9</th>\n",
       "      <td>liability</td>\n",
       "      <td>1</td>\n",
       "      <td>0</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "  product_type  rating  days\n",
       "0    financial       1     1\n",
       "1     property       1     0\n",
       "2    financial       1     0\n",
       "3    liability       1     0\n",
       "4     property       1     0\n",
       "5     property       0     1\n",
       "6    liability       1     0\n",
       "7     property       1     0\n",
       "8    financial       0     0\n",
       "9    liability       1     0"
      ]
     },
     "execution_count": 19,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "days = BayesianNodeRV('days', days_cpt, parent_names=['product_type', 'rating'])\n",
    "bn = FakeDataBayesianNetwork(product_type, rating)#, days)\n",
    "bn = FakeDataBayesianNetwork(product_type, rating, days)\n",
    "bn.rvs(10)"
   ]
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "b41a5a4101ae72c4a5f9743c63767ec11038be7e733d0952ee28fec8fab1bf9a"
  },
  "kernelspec": {
   "display_name": "Python 3.7.6 64-bit ('fake-data': virtualenv)",
   "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.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}