Truncated

Truncated.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import arviz as az\n",
    "import pymc3 as pm"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "3651\n"
     ]
    }
   ],
   "source": [
    "# create data\n",
    "np.random.seed(451)\n",
    "x = np.random.exponential(3, size=5000)\n",
    "\n",
    "minx=1\n",
    "maxx=20\n",
    "\n",
    "obs = x[np.where(~((x<minx) | (x>maxx)))] # remove values outside range\n",
    "print(len(obs))"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [λ]\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\nλ  2.967  0.051   2.872    3.063      0.001    0.001    3120.0  3118.0   \n\n   ess_bulk  ess_tail  r_hat  \nλ    3128.0    5648.0    1.0  ",
      "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>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_mean</th>\n      <th>ess_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>λ</th>\n      <td>2.967</td>\n      <td>0.051</td>\n      <td>2.872</td>\n      <td>3.063</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>3120.0</td>\n      <td>3118.0</td>\n      <td>3128.0</td>\n      <td>5648.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 3,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with pm.Model() as manual_model:\n",
    "    λ = pm.Exponential(\"λ\", lam=1/5)  # prior exponential with mean of 5\n",
    "    x = pm.Exponential('x', lam=1/λ, observed=obs) # obs exponential with mean of $\\lambda$.\n",
    "\n",
    "    exp_dist = pm.Exponential.dist(lam=1/λ) # this is not part of the model, just used to get the logcdf\n",
    "    norm_term = pm.Potential(\"norm_term\", -pm.math.logdiffexp(exp_dist.logcdf(maxx), exp_dist.logcdf(minx)) * x.size)\n",
    "\n",
    "    trace_manual= pm.sample(2000, tune=1000, return_inferencedata=True)\n",
    "\n",
    "az.summary(trace_manual)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [λ]\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 7 seconds.\n",
      "There were 40 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 46 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 40 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "There were 25 divergences after tuning. Increase `target_accept` or reparameterize.\n",
      "The rhat statistic is larger than 1.4 for some parameters. The sampler did not converge.\n",
      "The estimated number of effective samples is smaller than 200 for some parameters.\n"
     ]
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\nλ  2.733  0.432   2.109    3.365      0.216    0.165       4.0     4.0   \n\n   ess_bulk  ess_tail  r_hat  \nλ       4.0      11.0   3.56  ",
      "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>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_mean</th>\n      <th>ess_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>λ</th>\n      <td>2.733</td>\n      <td>0.432</td>\n      <td>2.109</td>\n      <td>3.365</td>\n      <td>0.216</td>\n      <td>0.165</td>\n      <td>4.0</td>\n      <td>4.0</td>\n      <td>4.0</td>\n      <td>11.0</td>\n      <td>3.56</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with pm.Model() as trunc_model:\n",
    "    λ = pm.Exponential(\"λ\", lam=1/5)\n",
    "    x = pm.Truncated(pm.Exponential, lower=minx, upper=maxx)('x', lam=1/λ, observed=obs)\n",
    "\n",
    "    trace_trunc = pm.sample(2000, tune=1000, return_inferencedata=True)\n",
    "az.summary(trace_trunc)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "outputs": [
    {
     "data": {
      "text/plain": "λ_log__      -1.00\nx         -7936.43\nName: Log-probability of test_point, dtype: float64"
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trunc_point = trunc_model.check_test_point(\n",
    "    test_point={'λ_log__': np.log(5)}\n",
    ")\n",
    "trunc_point"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "outputs": [
    {
     "data": {
      "text/plain": "λ_log__      -1.00\nx         -8749.23\nName: Log-probability of test_point, dtype: float64"
     },
     "execution_count": 6,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "manual_point = manual_model.check_test_point(\n",
    "    test_point={'λ_log__': np.log(5)}\n",
    ")\n",
    "manual_point"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "outputs": [
    {
     "data": {
      "text/plain": "λ_log__      0.0\nx         -812.8\nName: Log-probability of test_point, dtype: float64"
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "model_diff = manual_point - trunc_point\n",
    "model_diff"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "outputs": [
    {
     "data": {
      "text/plain": "array(812.80311981)"
     },
     "execution_count": 8,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "exp_dist = pm.Exponential.dist(1/5)\n",
    "norm_term = (-pm.math.logdiffexp(exp_dist.logcdf(maxx), exp_dist.logcdf(minx)) * x.size).eval()\n",
    "norm_term"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "outputs": [],
   "source": [
    "trunc_dist = pm.Truncated(pm.Exponential, lower=minx, upper=maxx).dist(lam=1/5)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "outputs": [
    {
     "data": {
      "text/plain": "-7936.430976274112"
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trunc_dist.logp(obs).eval().sum()"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "outputs": [
    {
     "data": {
      "text/plain": "-812.8031198123141"
     },
     "execution_count": 11,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "trunc_dist._normalization().eval() * len(obs)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [λ]\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 3 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\nλ  3.934  0.066   3.806    4.057      0.001    0.001    2535.0  2535.0   \n\n   ess_bulk  ess_tail  r_hat  \nλ    2533.0    5226.0    1.0  ",
      "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>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_mean</th>\n      <th>ess_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>λ</th>\n      <td>3.934</td>\n      <td>0.066</td>\n      <td>3.806</td>\n      <td>4.057</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>2535.0</td>\n      <td>2535.0</td>\n      <td>2533.0</td>\n      <td>5226.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 25,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with pm.Model() as trunc_model:\n",
    "    λ = pm.Exponential(\"λ\", lam=1/5)\n",
    "    x = pm.Truncated(pm.Exponential, lower=None, upper=None)('x', lam=1/λ, observed=obs)\n",
    "\n",
    "    trace_trunc = pm.sample(2000, tune=1000, return_inferencedata=True)\n",
    "az.summary(trace_trunc)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "outputs": [
    {
     "name": "stderr",
     "output_type": "stream",
     "text": [
      "Auto-assigning NUTS sampler...\n",
      "Initializing NUTS using jitter+adapt_diag...\n",
      "Multiprocess sampling (4 chains in 4 jobs)\n",
      "NUTS: [λ]\n",
      "Sampling 4 chains for 1_000 tune and 2_000 draw iterations (4_000 + 8_000 draws total) took 2 seconds.\n"
     ]
    },
    {
     "data": {
      "text/plain": "<IPython.core.display.HTML object>",
      "text/html": "\n    <div>\n        <style>\n            /* Turns off some styling */\n            progress {\n                /* gets rid of default border in Firefox and Opera. */\n                border: none;\n                /* Needs to be in here for Safari polyfill so background images work as expected. */\n                background-size: auto;\n            }\n            .progress-bar-interrupted, .progress-bar-interrupted::-webkit-progress-bar {\n                background: #F44336;\n            }\n        </style>\n      <progress value='0' class='' max='12000' style='width:300px; height:20px; vertical-align: middle;'></progress>\n      \n    </div>\n    "
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/plain": "    mean     sd  hdi_3%  hdi_97%  mcse_mean  mcse_sd  ess_mean  ess_sd  \\\nλ  3.937  0.064   3.819     4.06      0.001    0.001    3554.0  3554.0   \n\n   ess_bulk  ess_tail  r_hat  \nλ    3538.0    5514.0    1.0  ",
      "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>mean</th>\n      <th>sd</th>\n      <th>hdi_3%</th>\n      <th>hdi_97%</th>\n      <th>mcse_mean</th>\n      <th>mcse_sd</th>\n      <th>ess_mean</th>\n      <th>ess_sd</th>\n      <th>ess_bulk</th>\n      <th>ess_tail</th>\n      <th>r_hat</th>\n    </tr>\n  </thead>\n  <tbody>\n    <tr>\n      <th>λ</th>\n      <td>3.937</td>\n      <td>0.064</td>\n      <td>3.819</td>\n      <td>4.06</td>\n      <td>0.001</td>\n      <td>0.001</td>\n      <td>3554.0</td>\n      <td>3554.0</td>\n      <td>3538.0</td>\n      <td>5514.0</td>\n      <td>1.0</td>\n    </tr>\n  </tbody>\n</table>\n</div>"
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "with pm.Model() as manual_model:\n",
    "    λ = pm.Exponential(\"λ\", lam=1/5)  # prior exponential with mean of 5\n",
    "    x = pm.Exponential('x', lam=1/λ, observed=obs) # obs exponential with mean of $\\lambda$.\n",
    "    trace_manual= pm.sample(2000, tune=1000, return_inferencedata=True)\n",
    "\n",
    "az.summary(trace_manual)"
   ],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "outputs": [],
   "source": [],
   "metadata": {
    "collapsed": false,
    "pycharm": {
     "name": "#%%\n"
    }
   }
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 2
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython2",
   "version": "2.7.6"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 0
}