Added simple gaussian likelihood example

examples/simple.ipynb

@@ -0,0 +1,300 @@
+{
+ "cells": [
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "import numpy as np\n",
+    "import matplotlib.pylab as plt\n",
+    "import corner\n",
+    "import numpy as np\n",
+    "import glob\n",
+    "from PTMCMCSampler import PTMCMCSampler\n",
+    "\n",
+    "%matplotlib inline"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Define the likelihood and posterior\n",
+    "\n",
+    "Functions must read in parameter vector and output log-likelihood or log-prior. Usually easiest to use a class if you need to store some other data or parameters"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "class GaussianLikelihood(object):\n",
+    "    \n",
+    "    def __init__(self, ndim=2, pmin=-10, pmax=10):\n",
+    "        \n",
+    "        self.a = np.ones(ndim)*pmin\n",
+    "        self.b = np.ones(ndim)*pmax\n",
+    "        \n",
+    "        # get means\n",
+    "        self.mu = np.random.uniform(pmin, pmax, ndim)\n",
+    "\n",
+    "        # ... and a positive definite, non-trivial covariance matrix.\n",
+    "        cov  = 0.5-np.random.rand(ndim**2).reshape((ndim, ndim))\n",
+    "        cov  = np.triu(cov)\n",
+    "        cov += cov.T - np.diag(cov.diagonal())\n",
+    "        self.cov  = np.dot(cov,cov)\n",
+    "\n",
+    "        # Invert the covariance matrix first.\n",
+    "        self.icov = np.linalg.inv(self.cov)\n",
+    "        \n",
+    "    def lnlikefn(self, x):\n",
+    "        diff = x - self.mu\n",
+    "        return -np.dot(diff,np.dot(self.icov, diff))/2.0\n",
+    "    \n",
+    "    def lnpriorfn(self, x):\n",
+    "        \n",
+    "        if np.all(self.a <= x) and np.all(self.b >= x):\n",
+    "            return 0.0\n",
+    "        else:\n",
+    "            return -np.inf      \n"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup Gaussian model class"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "ndim = 20\n",
+    "pmin, pmax = 0.0, 10.0\n",
+    "glo = GaussianLikelihood(ndim=ndim, pmin=pmin, pmax=pmax)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Setup sampler\n",
+    "\n",
+    "Need to initalize the sample at ```p0``` and give an inital jump covariance matrix ```cov```."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# Set the start position and the covariance\n",
+    "p0 = np.random.uniform(pmin, pmax, ndim)\n",
+    "cov = np.eye(ndim) * 0.1**2"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "sampler = PTMCMCSampler.PTSampler(ndim, glo.lnlikefn, glo.lnpriorfn, np.copy(cov), outDir='./chains')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Add custom jump\n",
+    "\n",
+    "Can add custom jump in the following way"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "class UniformJump(object):\n",
+    "    \n",
+    "    def __init__(self, pmin, pmax):\n",
+    "        \"\"\"Draw random parameters from pmin, pmax\"\"\"\n",
+    "        self.pmin = pmin\n",
+    "        self.pmax = pmax\n",
+    "        \n",
+    "    def jump(self, x, it, beta):\n",
+    "        \"\"\" \n",
+    "        Function prototype must read in parameter vector x,\n",
+    "        sampler iteration number it, and inverse temperature beta\n",
+    "        \"\"\"\n",
+    "        \n",
+    "        # log of forward-backward jump probability\n",
+    "        lqxy = 0\n",
+    "        \n",
+    "        # uniformly drawm parameters\n",
+    "        q = np.random.uniform(self.pmin, self.pmax, len(x))\n",
+    "        \n",
+    "        return q, lqxy"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": [
+    "# add to jump proposal cycle\n",
+    "ujump = UniformJump(pmin, pmax)\n",
+    "sampler.addProposalToCycle(ujump.jump, 5)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Run Sampler for 100000 steps\n",
+    "\n",
+    "Different jump proposal weights are given as integers. For example we have used a weight of 20 for all three proposals here. That means that each will be used with a probability of 20/60 = 1/3."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false,
+    "scrolled": false
+   },
+   "outputs": [],
+   "source": [
+    "sampler.sample(p0, 100000, burn=500, thin=1, covUpdate=500,\n",
+    "               SCAMweight=20, AMweight=20, DEweight=20)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Get jump statistics\n",
+    "\n",
+    "Here you can track the acceptance rate of the different jump proposals used"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "jumpfiles = glob.glob('chains/*jump.txt')\n",
+    "jumps = map(np.loadtxt, jumpfiles)\n",
+    "for ct, j in enumerate(jumps):\n",
+    "    plt.plot(j, label=jumpfiles[ct].split('/')[-1].split('_jump.txt')[0])\n",
+    "plt.legend(loc='best', frameon=False)\n",
+    "plt.ylabel('Acceptance Rate')\n",
+    "plt.ylim(0.0, 1.1)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Get the data and plot the output\n",
+    "\n",
+    "The output data has ndim + 4 columns. The first ndim columns are just the samples from the parameters, the ndim+1 column is the log-posterior, ndim+2 is the log-likelihood, ndim+3 is the acceptance rate, and ndim+4 is the parallel tempering swap acceptance rate for the T=1 chain."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "data = np.loadtxt('chains/chain_1.txt')\n",
+    "chaint = data[:,:-4]"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# throw out first 25% of chain as burn in \n",
+    "burn = int(0.25*chaint.shape[0])\n",
+    "plt.plot(data[burn:,-4])"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "# get the true values and plot the posteriors for the first 10 parameters\n",
+    "truth = glo.mu\n",
+    "corner.corner(chaint[burn:,:10], bins=50, truths=truth);"
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": true
+   },
+   "outputs": [],
+   "source": []
+  }
+ ],
+ "metadata": {
+  "kernelspec": {
+   "display_name": "Python 2",
+   "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.10"
+  }
+ },
+ "nbformat": 4,
+ "nbformat_minor": 0
+}