From db5511ed663bee73d48d2c6d59f0523d215acdfd Mon Sep 17 00:00:00 2001 From: Justin Ellis Date: Tue, 23 Aug 2016 11:29:56 -0700 Subject: [PATCH] Added simple gaussian likelihood example --- examples/simple.ipynb | 300 ++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 300 insertions(+) create mode 100644 examples/simple.ipynb diff --git a/examples/simple.ipynb b/examples/simple.ipynb new file mode 100644 index 0000000..2748cc0 --- /dev/null +++ b/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 +}