C3/L9_numerical_stability.ipynb

{
  "nbformat": 4,
  "nbformat_minor": 0,
  "metadata": {
    "colab": {
      "name": "L8_classification_with_mixture_models.ipynb",
      "provenance": [],
      "collapsed_sections": [],
      "toc_visible": true
    },
    "kernelspec": {
      "name": "python3",
      "display_name": "Python 3",
      "language": "python"
    },
    "language_info": {
      "name": "python",
      "version": "3.8.5"
    }
  },
  "cells": [
    {
      "cell_type": "code",
      "metadata": {
        "id": "eDRfJeA9i540",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "caa4edea-6ba8-4638-cc88-876031272ee9"
      },
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt \n",
        "import pandas as pd \n",
        "import seaborn as sns\n",
        "import arviz as az\n",
        "from sklearn import datasets\n",
        "from statsmodels.tsa import stattools\n",
        "import statsmodels.api as sm\n",
        "import pymc3 as pm\n",
        "import pymc\n",
        "\n",
        "import sys\n",
        "sys.path.insert(0, '..')\n",
        "from utils.plot_lib import set_default\n",
        "\n",
        "set_default(figsize=(6, 4))"
      ],
      "execution_count": 1,
      "outputs": []
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "9c5CUB_VdYWL"
      },
      "source": [
        "### Numerical stability"
      ]
    },
    {
      "cell_type": "code",
      "metadata": {
        "id": "EV8NMNsvdQz0"
      },
      "source": [
        "from scipy.stats import norm\n",
        "\n",
        "## Consider a mixture of two normal distributions with equal weights (w1 = w2 = 1/2)\n",
        "## Component 1 has mean 0 and standard deviation 1\n",
        "## Component 2 has mean 1 and standard deviation 1\n",
        "## The observation is x = 50\n",
        "## What is Pr(c = 1 | x)?\n",
        "\n",
        "mu1 = 0\n",
        "mu2 = 1\n",
        "sig1 = 1\n",
        "sig2 = 1\n",
        "x = 50\n",
        "\n",
        "k1 = norm(loc = mu1, scale = sig1)\n",
        "k2 = norm(loc = mu2, scale = sig2)\n",
        "\n",
        "print('The probability that observation x = 50 belongs to component 1 is: {}'.format(k1.pdf(x)))\n",
        "print('The probability that observation x = 50 belongs to component 2 is: {}'.format(k2.pdf(x)))\n",
        "print('Computing Ratio. Numeral stability!: {}'.format(k1.pdf(x) / (k1.pdf(x) + k2.pdf(x))))"
      ],
      "execution_count": 2,
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "The probability that observation x = 50 belongs to component 1 is: 0.0\nThe probability that observation x = 50 belongs to component 2 is: 0.0\nComputing Ratio. Numeral stability!: nan\n"
          ]
        }
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 3,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.07585818002124355\n"
          ]
        }
      ],
      "source": [
        "## What if x=3?  Two ways to do the calculation\n",
        "## One way:  Direct calculation\n",
        "\n",
        "x = 3\n",
        "z1 = k1.pdf(3)\n",
        "z2 = k2.pdf(3)\n",
        "print(z1 / (z1+z2))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 4,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "0.07585818002124356\n"
          ]
        }
      ],
      "source": [
        "## A second way:  Compute in the logarithm scale, add b \n",
        "## to all values, and then exponentiate before standardizing\n",
        "\n",
        "## A second way:  Compute in the logarithm scale, add b \n",
        "## to all values, and then exponentiate before standardizing\n",
        "lz1 = k1.logpdf(x)\n",
        "lz2 = k2.logpdf(x)\n",
        "b = 3\n",
        "print(np.exp(lz1 - b) / (np.exp(lz1 - b) + np.exp(lz2 - b)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 5,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "3.1799709001977496e-22\n"
          ]
        }
      ],
      "source": [
        "# Other ways to compute it\n",
        "x = 50\n",
        "lz1 = k1.logpdf(x)\n",
        "lz2 = k2.logpdf(x)\n",
        "b = np.max([lz1, lz2])\n",
        "print(np.exp(lz1 - b) / (np.exp(lz1 - b) + np.exp(lz2 - b)))"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": 6,
      "metadata": {},
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "3.1799709001977496e-22\n"
          ]
        }
      ],
      "source": [
        "## Also right (just more cumbersome)\n",
        "lz1 = -0.5 * np.log(2 * np.pi) - 0.5 * x**2\n",
        "lz2 = -0.5 * np.log(2 * np.pi) - 0.5 * (x - mu2)**2\n",
        "b = np.max([lz1, lz2])\n",
        "print(np.exp(lz1 - b) / (np.exp(lz1 - b) + np.exp(lz2 - b)))"
      ]
    }
  ]
}