notebooks/example_python.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Example in Python\n",
    "\n",
    "This is a fairly minimal example, demonstrating the slightly different calling convention in the Python version of LRGS, compared with the R version.\n",
    "\n",
    "One notable and practical difference is that the Python version is significantly slower than the R version, despite using vectorization more efficiently in places. In the benchmarks I've done so far, it takes between 1.3 and 3 times longer to produce the same number of samples, depending on the model being fit.\n",
    "\n",
    "We'll just look at univariate regression, with p(x) a single Gaussian."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "import lrgs\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "One immediate difference is that we need to pass x and y as Numpy matrices, with each row corresponding to a data point. (So, in the univariate case, column vectors.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "## Use the same setup as the simple example in the R documentation\n",
    "## ... but with a bigger scatter and fewer data points, just because.\n",
    "x = np.matrix(np.random.normal(0.0, 5.0, size=100)).T\n",
    "y = np.matrix(np.random.normal(np.pi*x, np.exp(1.0)))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(x, y, 'o');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If we were going to pass measurement errors, it would be as a list of matrices, e.g."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# M = [np.asmatrix(np.eye(2)*1e-4) for i in range(x.shape[0])]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of doing everything in a single function call, like in R, we instead have one call to set up the model and a second to run the sampler. This does the first step (the first argument is the size of the Gaussian mixture):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "par = lrgs.Parameters_GaussMix(1, x, y)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "This object holds the current values of all the model parameters (currently rough default guesses) and has methods to update individual parameters or all parameters. We could call these directly to sample the parameters \"in place\", but of course we actually want to store the chain. For that, we use a second declaration. Here the second argument is the length of the chain we want to run, and the third encodes the parameters that we want to store (in the usual LRGS way):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "chain = lrgs.Chain(par, 100, 'bsmt')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In the \"run\" call, we specify any parameters to be fixed (x and y in this case):"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "chain.run(fix='xy')"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The chain can now be converted to a dict or recarray (and further converted to a pandas data frame, for e.g.). Note that the \"proper names\" of the parameters to include must be specified here, and that the resulting keys differ slightly from the column naming scheme used in the R function `Gibbs.post2dataframe`. (In particular, indices are offset by underscores and start from zero, consistent with Python indexing.)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "dchain = chain.to_dict([\"B\", \"Sigma\", \"mu\", \"Tau\"])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, `dchain[\"B_1_0\"]` is the slope; in R it would have been called `B21`."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Plot things. Red lines show the input values."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "keep = range(10, chain.nmc) # throw out a little burn-in"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"B_0_0\"][keep], 'o');\n",
    "plt.axhline(y=0.0, color='r');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"B_1_0\"][keep], 'o');\n",
    "plt.axhline(y=np.pi, color='r');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"B_0_0\"][keep], dchain[\"B_1_0\"][keep], 'o');\n",
    "plt.plot(0.0, np.pi, 'or');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"Sigma_0_0\"][keep], 'o');\n",
    "plt.axhline(y=np.exp(2.0), color='r');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"mu_0_0\"][keep], 'o');\n",
    "plt.axhline(y=0.0, color='r');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "plt.plot(dchain[\"Tau_0_0_0\"][keep], 'o');\n",
    "plt.axhline(y=5.0**2, color='r');"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python [default]",
   "language": "python",
   "name": "python2"
  },
  "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.12"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}