Ch2.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Notes\n",
    "---\n",
    "\n",
    "## Equations:\n",
    "\n",
    "### Gausian (Normal) Distribution\n",
    "- $\\mathcal{N}(x;\\mu,\\sigma^2)$\n",
    "- $p(x) = (2\\pi\\sigma^{2})^{-\\frac{1}{2}}\\ exp\\{-\\frac{1}{2} \\frac{(x-\\mu)^2}{\\sigma^2} \\}$\n",
    "- $p(x) = det(2\\pi\\Sigma)^{-\\frac{1}{2}}\\ exp\\{-\\frac{1}{2} (x-\\mu)^T \\Sigma^-1 (x-\\mu) \\}$ (Multivariate Distribution)\n",
    "\n",
    "### Independence\n",
    "- $p(x,y) = p(x)\\ p(y)$\n",
    "- $p(x\\ |\\ y) = \\frac{p(x,y)}{p(y)} = \\frac{p(x)\\ p(y)}{p(y)} = p(x)$\n",
    "\n",
    "### Conditional Independence\n",
    "- When y does not cary information about x if z is known:\n",
    "- $p(x\\ |\\ z) = p(x\\ |\\ z,y)$\n",
    "- $p(y\\ |\\ z) = p(y|z,x)$\n",
    "- This does not imply absolute indepence of x and y\n",
    "\n",
    "### Total Probability\n",
    "- $p(x) = \\Sigma_y p(x\\ |\\ y)p(y)$ (discrete)\n",
    "- $p(x) = \\int p(x\\ |\\ y)p(y) dy$ (continuous)\n",
    "\n",
    "### Bays Rule\n",
    "- $p(y\\ |\\ x) = \\frac{p(x\\ |\\ y)p(x)}{p(y)} = \\frac{p(x\\ |\\ y)p(x)}{\\Sigma_{x'}p(y\\ |\\ x')p(x')}$ (discrete)\n",
    "- $p(y\\ |\\ x) = \\frac{p(x\\ |\\ y)p(x)}{p(y)} = \\frac{p(x\\ |\\ y)p(x)}{\\int p(y\\ |\\ x')p(x') dx'}$ (continuous)\n",
    "- $p(y\\ |\\ x) = \\eta \\space p(x\\ |\\ y)p(x)$ ($\\eta$ is a normalizer)\n",
    "- $p(x\\ |\\ y,z) = \\frac{p(y\\ |\\ x,z)p(x\\ |\\ z)}{p(y\\ |\\ x)}$\n",
    "- $p(x,y\\ |\\ z) = p(x\\ |\\ z)p(y\\ |\\ z)$\n",
    "\n",
    "### Expectation\n",
    "- $E[X] = \\Sigma_x x \\space p(x)$ (discrete)\n",
    "- $E[X] = \\int x \\space p(x) \\space dx$ (continuous)\n",
    "- $E[aX+b] = a \\space E[X] + b$\n",
    "\n",
    "### Covariance\n",
    "- The squared expected deviatino from the mean\n",
    "- $Cov[X] = E[X - E[X]]^2 = E[X^2] - E[X]^2$\n",
    "\n",
    "### Entropy\n",
    "- The expected informationh that the value x carries\n",
    "- $H_p (x) = E[-log_2(p(x))]$\n",
    "- $H_p (x) = -\\Sigma_x p(x) log_2(p(x))$\n",
    "- $H_p (x) = -\\int p(x) \\space log_2(p(x)) \\space dx$\n",
    "\n",
    "---\n",
    "## Acronyms\n",
    "- PDF = Probability Density Function\n",
    "- HMM = Hidden Markov Model\n",
    "- DBN = Dynamic Bayes Network\n",
    "\n",
    "### Definitions\n",
    "\n",
    "---\n",
    "## Robot Environment Interaction\n",
    "\n",
    "- State = $x$ - Collection of all aspects of the robot and the environment that can affect the future\n",
    "    - $x_t$ is the state $x$ at time $t$\n",
    "- Pose - Robot location, orientation, and configuration\n",
    "- Control Action = $u$\n",
    "    - Control Actions decreases information about the state\n",
    "- Measurement = $z$\n",
    "    - Measurement increases information about the state\n",
    "- Belief = $bel(x_t)$ - Belief about the state $x$ at time $t$\n",
    "    - $\\overline{bel}(x_t)$ - Belief about the state $x$ at time $t$ before incorperating measurement $z_t$\n",
    "- State Update\n",
    "    - Assume state $x_t$ provides sufficient information about all previous states, control actions, and measurements\n",
    "    - Assume the robot makes a control action, then takes a measurement\n",
    "    - State Transition Probability $p(x_t\\ |\\ x_{t-1}, u_t)$\n",
    "    - Measurement Probability $p(z_t\\ |\\ x_{t})$\n",
    "- Probability Distribuion\n",
    "    - Prediction $\\overline{bel}(x_t)$\n",
    "        - Predicts the state after Control Action\n",
    "    - Correction $bel(x_t)$\n",
    "        - Corrects state estimation using measurement information\n",
    "\n",
    "---\n",
    "\n",
    "## Bayes Filter\n",
    "\n",
    "1. **Algorithm Bayes_filter**($bel(x_t-1), u_t, z_t$):\n",
    "1. $\\ \\ $for all $x_t$ do\n",
    "1. $\\ \\ $$\\ \\ $$\\overline{bel}(x_t) = \\int\\ p(x_t\\ |\\ u_t, x_{t-1})\\ bel(x_{t-1})\\ dx_{t-1}$\n",
    "1. $\\ \\ $$\\ \\ $$bel(x_t) = \\eta\\ p(z_t\\ |\\ x_t)\\ \\overline{bel}(x_t)$\n",
    "1. $\\ \\ $end for\n",
    "1. $\\ \\ $return $bel(x_t)$\n",
    "\n",
    "\n",
    "- Step 3 is prediction\n",
    "- Step 4 is the measurement update (correction)\n",
    "- The value of $bel(x_t)$ is fed back into the algorithm for the next iteration\n",
    "- The final result is normaized to 1 by normalization constant $\\eta$\n",
    "- $bel(x_0)$ must be given\n",
    "    - If is is known, it should be a point mass\n",
    "    - If it is NOT known, it should be a uniform distribution\n",
    "    - If partial knowledge, it should be a some non-uniform distribution (not common in practice)\n",
    "- This simple form only works for very simple estimation problems\n",
    "\n",
    "---\n",
    "## Markov Assumption\n",
    "- The past and future states are independent if you know the current state $x_t$\n",
    "- Will not be true for almost any problem we work on, but the approximation is generally good enough"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercises"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 1. Robot Range Finder\n",
    "\n",
    "- Range finder measures between 0m and 3m\n",
    "- Assume the probability of measureing distances is evently distributed\n",
    "- If the range finder is faultly, it will always measures < 1m\n",
    "- Prior probability of sensor being faultly is p=0.01\n",
    "- If the sensor is meausered N times and the value is always below 1m, what is the posterior probability that the sensor is faulty?\n",
    "- Use Bayes's theorm:\n",
    "    - $p(bad\\_sensor\\ |\\ measurement) = \\frac{p(measurement\\ |\\ bad\\_sensor)\\ p(bad\\_sensor)}{p(measurement)}$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr style=\"text-align: right;\">\n",
       "      <th>Samples</th>\n",
       "      <th>0</th>\n",
       "      <th>1</th>\n",
       "      <th>2</th>\n",
       "      <th>3</th>\n",
       "      <th>4</th>\n",
       "      <th>5</th>\n",
       "      <th>6</th>\n",
       "      <th>7</th>\n",
       "      <th>8</th>\n",
       "      <th>9</th>\n",
       "      <th>10</th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th>p(bad|N measurements)</th>\n",
       "      <td>0.0100</td>\n",
       "      <td>0.0294</td>\n",
       "      <td>0.0833</td>\n",
       "      <td>0.2143</td>\n",
       "      <td>0.4500</td>\n",
       "      <td>0.7105</td>\n",
       "      <td>0.8804</td>\n",
       "      <td>0.9567</td>\n",
       "      <td>0.9851</td>\n",
       "      <td>0.9950</td>\n",
       "      <td>0.9983</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "Samples                   0      1      2      3      4      5      6      7   \\\n",
       "p(bad|N measurements) 0.0100 0.0294 0.0833 0.2143 0.4500 0.7105 0.8804 0.9567   \n",
       "\n",
       "Samples                   8      9      10  \n",
       "p(bad|N measurements) 0.9851 0.9950 0.9983  "
      ]
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "import numpy as np\n",
    "import pandas as pd\n",
    "pd.set_option('float_format', lambda x: '{:.4f}'.format(x))\n",
    "\n",
    "# Initialize the constants\n",
    "prior_prob_bad = 0.01\n",
    "prob_meas_lt_1_given_bad = 1\n",
    "prob_meas_lt_1_given_good = 1 / 3\n",
    "\n",
    "# Initialize the probability that the sensor is bad\n",
    "prob_bad = prior_prob_bad\n",
    "\n",
    "probabilites = []\n",
    "\n",
    "# Loop over the measurements\n",
    "for N in range(10+1):\n",
    "    \n",
    "    # Save the probability\n",
    "    probabilites.append(prob_bad)\n",
    "    \n",
    "    # Compute the probability of the measurement being less than 1\n",
    "    prob_good = 1 - prob_bad  # Probablity the sensor is good\n",
    "    prob_meas = prob_bad  * prob_meas_lt_1_given_bad + \\\n",
    "                prob_good * prob_meas_lt_1_given_good\n",
    "\n",
    "    # Update the probability that the sensor is bad\n",
    "    prob_bad = (prob_meas_lt_1_given_bad * prob_bad) / prob_meas\n",
    "\n",
    "\n",
    "# Print the result\n",
    "result = pd.DataFrame(data=np.array([probabilites]))\n",
    "result.index = ['p(bad|N measurements)',]\n",
    "result.columns.name = 'Samples'\n",
    "result"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 2. Weather simulation\n",
    "\n",
    "- The weather can take on three values, sunny, cloudy, and rainy\n",
    "- The weather transition function is a Markov chain with this transition table\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<div>\n",
       "<table border=\"1\" class=\"dataframe\">\n",
       "  <thead>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>when</th>\n",
       "      <th colspan=\"3\" halign=\"left\">tomorrow</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th></th>\n",
       "      <th>weather</th>\n",
       "      <th>sunny</th>\n",
       "      <th>cloudy</th>\n",
       "      <th>rainy</th>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>when</th>\n",
       "      <th>weather</th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "      <th></th>\n",
       "    </tr>\n",
       "  </thead>\n",
       "  <tbody>\n",
       "    <tr>\n",
       "      <th rowspan=\"3\" valign=\"top\">today</th>\n",
       "      <th>sunny</th>\n",
       "      <td>0.8000</td>\n",
       "      <td>0.2000</td>\n",
       "      <td>0.0000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>cloudy</th>\n",
       "      <td>0.4000</td>\n",
       "      <td>0.4000</td>\n",
       "      <td>0.2000</td>\n",
       "    </tr>\n",
       "    <tr>\n",
       "      <th>rainy</th>\n",
       "      <td>0.2000</td>\n",
       "      <td>0.6000</td>\n",
       "      <td>0.2000</td>\n",
       "    </tr>\n",
       "  </tbody>\n",
       "</table>\n",
       "</div>"
      ],
      "text/plain": [
       "when          tomorrow              \n",
       "weather          sunny cloudy  rainy\n",
       "when  weather                       \n",
       "today sunny     0.8000 0.2000 0.0000\n",
       "      cloudy    0.4000 0.4000 0.2000\n",
       "      rainy     0.2000 0.6000 0.2000"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "transition_table = pd.DataFrame(np.array([[.8, .2, 0], [.4, .4, .2], [.2, .6, .2]]))\n",
    "transition_table.columns=pd.MultiIndex.from_arrays([['tomorrow']*3, ['sunny', 'cloudy', 'rainy']])\n",
    "transition_table.columns.names = ['when', 'weather']\n",
    "transition_table.index = pd.MultiIndex.from_arrays([['today']*3, ['sunny', 'cloudy', 'rainy']])\n",
    "transition_table.index.names = ['when', 'weather']\n",
    "transition_table"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (a)\n",
    "- Suppose day 1 = sunny\n",
    "- Find the probability of the following sequnce of days:\n",
    "    - Day 2 = cloudy\n",
    "    - Day 3 = cloudy\n",
    "    - Day 4 = rainy"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "probablity of sequence = 1.60%\n"
     ]
    }
   ],
   "source": [
    "todays_weather = 'sunny'\n",
    "forcast = ['cloudy', 'cloudy', 'rainy']\n",
    "prob_of_sequence = 1\n",
    "\n",
    "for tomorrows_weather in forcast:\n",
    "    prob_of_transition = transition_table['tomorrow', tomorrows_weather].loc['today', todays_weather]\n",
    "    prob_of_sequence *= prob_of_transition\n",
    "    todays_weather = tomorrows_weather\n",
    "\n",
    "print('probablity of sequence = {:.2f}%'.format(prob_of_sequence*100))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (b)\n",
    "- Write a simulator that can randomly generate sequences of \"weather\" using this state transition matrix"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "['cloudy', 'cloudy', 'sunny', 'cloudy', 'cloudy', 'sunny', 'sunny', 'sunny', 'sunny', 'sunny']\n"
     ]
    }
   ],
   "source": [
    "import random\n",
    "\n",
    "class WeatherSimulator:\n",
    "    '''Class that simulates the weather given a transition table and starting weather.'''\n",
    "    \n",
    "    def __init__(self, transition_table, start_weather):\n",
    "        self.todays_weather = start_weather\n",
    "        \n",
    "        # Precomute the cumulative transition probability\n",
    "        cum_probs = dict()\n",
    "        for weather in transition_table.loc['today'].index: # all possible weather states\n",
    "            next_weather = transition_table.loc['today'].loc[weather]['tomorrow'] # all possible transitions\n",
    "            next_weather_states = next_weather.index.values # names (eg. sunny)\n",
    "            next_weather_cum_probs = np.cumsum(next_weather).values # cumulative sum of probabilities\n",
    "            cum_probs[weather] = (next_weather_states, next_weather_cum_probs)\n",
    "        self.cum_probs = cum_probs\n",
    "        \n",
    "    def simulate(self):\n",
    "        '''Choses a next weather state at random given the transition probabilities.\n",
    "            The simulator state is updated to this new state and this new state is returned.'''\n",
    "        possible_weather, cumulative_sum = self.cum_probs[self.todays_weather]\n",
    "        r = random.random()\n",
    "        #print(r, cumulative_sum, cumulative_sum > r, possible_weather[cumulative_sum > r])\n",
    "        tomorrows_weather = possible_weather[cumulative_sum > r][0]\n",
    "        self.todays_weather = tomorrows_weather\n",
    "        return tomorrows_weather\n",
    "    \n",
    "ws = WeatherSimulator(transition_table, 'sunny')\n",
    "print([ws.simulate() for i in range(10)])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (c)\n",
    "- Use the simulator to come up with the steady state weather probabilities"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p(sunny) = 64.26%\n",
      "p(cloudy) = 28.64%\n",
      "p(rainy) = 7.10%\n"
     ]
    }
   ],
   "source": [
    "# sunny, cloudy, rainy\n",
    "weather_states = transition_table.loc['today'].index.values\n",
    "\n",
    "# Initalize the results dict\n",
    "results = dict()\n",
    "for weather_state in weather_states:\n",
    "    results[weather_state] = 0\n",
    "\n",
    "# Initalize the weather simulator\n",
    "ws = WeatherSimulator(transition_table, 'sunny')\n",
    "\n",
    "# Run N iterations and count the times the weather was in each state\n",
    "N = 100000\n",
    "for i in range(N):\n",
    "    results[ws.simulate()] +=1\n",
    "\n",
    "# Print the results\n",
    "for weather in weather_states:\n",
    "    print('p({}) = {:.2f}%'.format(weather, results[weather]*100.0/N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (d)\n",
    "- Find the closed form solution for the steady state weather probabilities\n",
    "- We can write the system of equations:\n",
    "    - $p(sunny) = \\Sigma_{weather}\\ p(sunny\\ |\\ weather)$\n",
    "        - S = (.8S + .4C + .2R)\n",
    "        - 0 = -.2S + .4C + .2 R\n",
    "    - $p(cloudy) = \\Sigma_{weather}\\ p(cloudy\\ |\\ weather)$\n",
    "        - C = (.2S + .4C + .6R)\n",
    "        - 0 = .2S + -.6C + .6R\n",
    "    - 1 = $\\Sigma_{weather}\\ p(weather)$\n",
    "        - 1.0 = S + C + R\n",
    "- Sovling this systems of equations for S,C,R will yeild the steady state probabilities\n",
    "        "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "p(sunny) = 64.26%\n",
      "p(cloudy) = 28.64%\n",
      "p(rainy) = 7.10%\n"
     ]
    }
   ],
   "source": [
    "weather_states = ['sunny', 'cloudy', 'rainy']\n",
    "\n",
    "a = np.array([[-0.2, 0.4, 0.2],\n",
    "              [0.2, -0.6, 0.6],\n",
    "              [1.0, 1.0, 1.0]])\n",
    "b = np.array([0.0,0.0,1.0])\n",
    "\n",
    "# Solve the system of equations\n",
    "p_sunny, p_cloudy, p_rainy = np.linalg.solve(a,b)\n",
    "\n",
    "# Print the results\n",
    "for weather in weather_states:\n",
    "    print('p({}) = {:.2f}%'.format(weather, results[weather]*100.0/N))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### (e)\n",
    "- Solve for the entropy of this system\n",
    "- entropy = -∑p(x)*log2(p(x))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "collapsed": false
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "entropy = 1.1981\n"
     ]
    }
   ],
   "source": [
    "# Grab the results from (d)\n",
    "probabilities = np.array([p_sunny, p_cloudy, p_rainy])\n",
    "\n",
    "entropy = -sum(probabilities * np.log2(probabilities))\n",
    "\n",
    "print('entropy = {:.4f}'.format(entropy))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 3. Noisy Sensor"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# [sunny, cloudy, rainy]\n",
    "p_w_0 = np.array([1.0,0.0,0.0])\n",
    "\n",
    "def weather_update(w):\n",
    "    \"\"\"Returns the probabiliy of tomorrows weather given the probablity to todays weather.\"\"\"\n",
    "    t = np.array([[0.8, 0.2, 0.0],\n",
    "                  [0.4, 0.4, 0.2],\n",
    "                  [0.2, 0.6, 0.2]])\n",
    "    return np.dot(w,t)\n",
    "\n",
    "def normalize(p):\n",
    "    \"\"\"Normailzes the probablities to they add up to One\"\"\"\n",
    "    return p / p.sum()\n",
    "\n",
    "def measure_update(w, m):\n",
    "    \"\"\"Returns the probabiliy of todays weather given the prior probablity to todays weather and a measurement.\"\"\"\n",
    "    t = np.array([[0.6, 0.4, 0.0],\n",
    "                  [0.3, 0.7, 0.0],\n",
    "                  [0.0, 0.0, 1.0]]).T\n",
    "    \n",
    "    raw = w*t[m]\n",
    "    return raw / raw.sum()\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ms = [1,1,2,0]\n",
    "w = p_w_0[:]\n",
    "for i, m in enumerate(ms):\n",
    "    print 'day %d, m = %d' % (i + 2, m)\n",
    "    w = weather_update(w)\n",
    "    print '\\t prior = %s' % str(w)\n",
    "    w = measure_update(w,m)\n",
    "    print '\\t post = %s' % str(w)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 4. Continuous Estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "def gauss(mu, sigma_sq, x):\n",
    "    return (2*np.pi*sigma_sq)**-0.5 *np.exp(-(x-mu)**2/(2*sigma_sq))\n",
    "\n",
    "gauss(0,np.sqrt(0.1),0.1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "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.5.1"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}