day5_exercises.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "# Day 5: Integrated auto-correlation time and Bayesian posterior sampling\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For today's tutorial you will need to import the following modules"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "# import numpy and matplotlib\n",
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "import scipy.stats as scstats\n",
    "# import mini MD modules\n",
    "from miniMD import models\n",
    "from miniMD import integrators\n",
    "from miniMD import outputshedulers as outp\n",
    "# import function to compute autocorrelation\n",
    "from miniMD.integrators import autocorr"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Exercise 1\n",
    "The code below creates $N=100000$ independent samples drawn from a standard normal distribution $\\mu(q) = \\rho(q) d q$,\n",
    "$$\n",
    "\\rho(q) = \\frac{1}{z}e^{-U(q)}\n",
    "$$\n",
    "with $U(q) =\\frac{1}{2}q^2$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Create iid samples  \n",
    "Nsteps = 100000\n",
    "xiid = np.random.normal(0., 1., Nsteps)\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The following code snippet bins the samples generate above, and plots the corresponding histogram against the exact PDF of a normal distsribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# Bin samples\n",
    "qbins = np.linspace(-3,3,50)\n",
    "Hq = np.histogram(xiid,bins=qbins,density=True)[0]\n",
    "deltaq = qbins[1]-qbins[0]\n",
    "qcenter = .5*deltaq + qbins[:-1]\n",
    "\n",
    "#Plot results\n",
    "fig, ax = plt.subplots()\n",
    "ax.plot(qcenter,Hq,label=\"Empirical\")\n",
    "ax.plot(qcenter, scstats.norm.pdf(qcenter),label=\"Exact\")\n",
    "ax.legend()\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Repeat the above experiment but now use a samples generated using the BAOAB Langevin sampler (you can use the code below). What do you observe?\n",
    "- Create trace plots of the trajectories both for the iid sample and the sample generated with BAOAB. What do you observe?\n",
    "- Change the values of the friction_constant and the stepsize. How does this effect appearance of the trace plot and the empirical distribution?\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "# initialise the model\n",
    "model = models.HarmonicOscillator()\n",
    "# Sample using BAOAB\n",
    "integrator = integrators.LangevinBAOAB(model=model, \n",
    "                         h=.1, \n",
    "                         Tk_B=1.0,\n",
    "                         gamma=10, \n",
    "                        ) \n",
    "op = outp.BufferedOutputsheduler(integrator, \n",
    "                                 Nsteps=20000, \n",
    "                                 varname_list=['q','p'], \n",
    "                                 modprnt=1)\n",
    "integrator.run(initial_values={'q' : np.zeros(model.dim),'p' : np.zeros(model.dim)}) # sample the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You should have observed in the above exercise that (at least for some parametrizations of the BAOAB method) the histogram generated with the BAOAB method is much more wiggly than the histogram generated using i.i.d. samples. The reason for that is that samples generated with the BAOAB integrator are correlated, which results in a reduced \"effective\" sample size, thus more \"wigglyness\" in the histogram. This observation is made more precise by relating the so called integrated autocorrelation time of the process to the variance in the central limit theorem for the Monte-Carlo estimate:\n",
    "\n",
    "### Central limit Theorem for underdamed Langevin dynamics\n",
    "Let\n",
    "$$\n",
    "\\mathbb{E}_{\\rho}\\varphi := \\int \\varphi(q) \\rho(q) d q\n",
    "$$\n",
    "and \n",
    "$$\n",
    "\\sigma_\\varphi^2 = \\int\\left ( \\varphi(q)-\\mathbb{E}_{\\rho}\\varphi \\right )^2 \\rho(q) dq\n",
    "$$\n",
    "denote the expectation and variance of the observable $\\varphi$ with respect to the probability density $\\rho$, respectively.\n",
    "\n",
    "Under mild conditions on the potential function $U$, one can show that for large $T>0$ a central limit theorem of the form \n",
    "\n",
    "$$\n",
    "\\sqrt{T}\\left ( \\frac{1}{T}\\int_0^T \\varphi(q(t)) d t - \\mathbb{E}_{\\rho}\\varphi  \\right )  \\sim \\mathcal{N}(0,\\sigma_\\varphi^2 \\tau_\\varphi)\n",
    "$$\n",
    "\n",
    "holds for the stochastic dynamics of the underdamped Langevin equation, where $\\tau_\\varphi$ denotes the so called integrated autocorrelation time of the stochastic process $\\varphi(q(t))$, i.e.,\n",
    "\n",
    "$$\n",
    "\\tau_\\varphi = \\int_0^{\\infty}C_{\\varphi}(t) d t\n",
    "$$\n",
    "\n",
    "where\n",
    "\n",
    "$$\n",
    "C_{\\varphi}(t) = \\mathbb{E}\\left [\\left (\\varphi(q(t))-\\mathbb{E}_{\\rho}\\varphi \\right )\\left (\\varphi(q(0))-\\mathbb{E}_{\\rho}\\varphi \\right )  \\right ], q(0) \\sim \\rho.\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Compute for the above Harmonic Oscillator mode the autocorrelation function for the observable $\\varphi(q) = q^2$ and trace plots both for $q$ and $q^2$. You can use the below code snipped to compute the autocorrelation for this obervable up to a maximum lag time specified by the argument maxlag. Can you characterise how the form of the autocorrelation function relates to the statistical properties of the sample?\n",
    "- For fixed stepsize $h$, can identify the optimal choice for the parameter gamma leading to the smallest integrated autocorrelation time."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "fig,ax = plt.subplots()\n",
    "autocfunc = autocorr(op.traj_q**2,maxlag=100)\n",
    "ax.plot(autocfunc,label=\"Autocorrelation function\")\n",
    "ax.axhline(y=0, color='k')\n",
    "ax.axvline(x=0, color='k')\n",
    "ax.legend()\n",
    "plt.show()\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 2\n",
    "In this exercise we will explore how the Ensemble Quasi Newton scheme which was introduced in today's lecture (see also Matthews, Weare, Leimkuhler - 2016 - Ensemble preconditioning for Markov chain Monte Carlo simulation or\n",
    "https://link.springer.com/article/10.1007/s11222-017-9730-1) can improve the sampling in the case of a badly conditioned target density. We will focus on the case of a Gaussian target density, i.e.,  \n",
    "\n",
    "$$\n",
    "\\rho_{target}(x) \\propto e^{- \\frac{1}{2}x^{T}\\Omega x}\n",
    "$$\n",
    "\n",
    "where the matrix $\\Omega$ has a large condition number,i.e., the ratio between the largest eigenvalue of $\\Omega$ and the smallest eigenvalue of $\\Omega$ is \"large\". For $\\Omega$ being a diagonal matrix, such a target density can be specified using the code below:\n",
    "\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "#Specify precision matrix of multivariate Gaussian target density \n",
    "M = 1000.0\n",
    "k_vec = np.array([1.0,1/M])\n",
    "model = models.MVHarmonicOscillator(k_vec)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "As you learned in the lecture, the Ensemble Quasi Newton method operates on an ensemble of replicas of the system under consideration. In order to create such an ensemble of replicas for the model you specified above you can make use of the class ReplicatedModel which is implimented in './miniMD/models.py' as a subclass of the class Model:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "replicated_model = models.ReplicatedModel(model, # model specifying target density \n",
    "                                          nreplicas=10 # number of replicas/copies to be created\n",
    "                                         )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A version of the Ensemble Quasi Newton method (using global covariance estimates) is implemented as EnsembleQuasiNewton as a subclass of the class LangevinBAOSplitting within the module miniMD.integrators. This sampler/integrator can be initialized by providing an instance of the class ReplicatedModel and the usual parameters for Langevin samplers as well as the parameter `B_update_mod`, which specifies the lag time between updates of the preconidioner materix $B$:"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "eqn_sampler = integrators.EnsembleQuasiNewton(repmodel=replicated_model, \n",
    "                            h=.1, \n",
    "                            Tk_B=1.0,\n",
    "                            gamma=1.0, \n",
    "                            B_update_mod=10\n",
    "                            ) "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After initialising \"eqn_sampler\" you can use this object to sample from your target density, e.g.,"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "op = outp.BufferedOutputsheduler(eqn_sampler, \n",
    "                                 Nsteps=10000, \n",
    "                                 varname_list=['q','p'], \n",
    "                                 modprnt=1)\n",
    "eqn_sampler.run(initial_values={'q' : np.random.normal(0.0,1.0,replicated_model.dim),'p' : np.zeros(replicated_model.dim)}) # sample the system"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For fixed stepsize $\\Delta t = .1$ compare the performance of a standard Langevin BAOAB integrator (Hint: for nreplicas=1, the EnsembleQuasiNewton method is identical to standard BAOAB) and the Ensemble Quasi Newton method for the 2-dimensional Multivariate Gaussian specified above by inspecting trace plots of $x_1$ and $x_2$ and/or the autocorrelation function of these variables.\n",
    "-  vary the condition number of $\\Omega$ by changing the value of $M$ (Hint: better only increase M, so that you don't run into stability issues). How does this affect the performance of standard BAOAB and the Ensemble Quasi Newton method, respectively?\n",
    "- How sensitive is the performance of the Ensemble Quasi Newton method on the number of replicas? \n",
    "- __Optional__: You can experiment with changing the frequency at which the preconditioning matrix is updated by modifying the optional argument \"B_update_mod\" when initializing the EQN sampler (The deault value is B_update_mod=1). "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Exercise 3"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "Consider a binary classification problem, for a given set of paired data $(x_i,y_i)_{1 \\leq i \\leq N}$, where $y_i \\in \\mathbb{R}^d$ denotes the predictor variable and $y_i \\in \\{0, 1\\}$ the response variable. Let $X_i$ and $Y_i$ as described above be i.i.d. random variables associated with the observations $(y_i,z_i)$, respectively. A common way to model the dependence between the predictor variable and the response variable in binary classification problems is to assume that the conditional probability of the event $Y_i = 1$ given $X_i$ is described by the values of the logistic function evaluated at a linear transformation of $X_i$, i.e.,\n",
    "\n",
    "$$\n",
    "\\mathbb{P}(Y_i =1|X_i)=f( x_i \\beta^T\\,)\n",
    "$$\n",
    "\n",
    "hence \n",
    "$$\\mathbb{P}(Y_i \\beta^T X_i) = f(Y_i \\beta^T X_i),$$\n",
    "where\n",
    "$$\n",
    "f(s) = \\frac{1}{1+e^{-s}}.\n",
    "$$\n",
    "is the logistic function and the weights w ∈ Rn define the linear transformation. In a Bayesian setup, the weights are considered to be random variables distributed according to a prior distribution, w ∼ π. By Bayes rule, the posterior distribution over w is then found to be,\n",
    "$$\n",
    "p(\\beta | x,y) \\propto \\pi(\\beta) p(y, | x, \\beta) = \\pi(\\beta) \\prod_{i=1}^N f( x_i \\beta^T\\, y_i ). \n",
    "$$\n",
    "Using the posterior distribution we can predict the class label of a new data point as \n",
    "$$\n",
    "\\mathbb{P}[y_{pred} | x_{pred},x,y] = \\int f(x_{pred}^T\\beta) p(\\beta | x,y) d \\beta\n",
    "$$\n",
    "In practice the integral on the left hand side of the equation is intractable. This is were sampling comes into play. One generate a sample $(\\beta^(k))_{1\\leq k \\leq N_{steps}}$ from the posterior distribution and approximate the integral term as \n",
    "$$\n",
    "\\int f(x_{pred}^T\\beta) p(\\beta | x,y) d \\beta  \\approx \\frac{1}{N_{steps}}\\sum_{k=1}^{N_{steps}} f(x_{pred}^T\\beta^{(k)})\n",
    "$$\n",
    "\n",
    "In this exercise we will consider a simple logistic regression classfier using a synthetic dataset of $N$ paired observations $(x_i,y_i)_{1\\leq i \\leq N}$, with $x_i \\in \\mathbb{R}^2$ being the predictor variable and $y_i\\in \\{0,1\\}$ being the class label as a training set. The code below creates such a synthetic dataset."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": [
    "np.random.seed(seed=11) #Fix seed \n",
    "\n",
    "data_dim = 2 # Dimension of predictor variable\n",
    "Ndata1 = 10  # Number of points with class label 0\n",
    "Ndata2 = 10  # Number of points with class label 1\n",
    "\n",
    "mu1 = np.array([-4,0]) # mean of predictor variables with class label 0\n",
    "mu2 =  np.array([4,0]) # mean of predictor variables with class label 1\n",
    "cov1 = np.eye(data_dim) # covariance of predictor variables with class label 0\n",
    "cov2 = np.eye(data_dim) # covariance of predictor variables with class label 1\n",
    "\n",
    "\n",
    "# Sample data points \n",
    "X1 = np.random.multivariate_normal(mu1,cov1,size=Ndata1)\n",
    "Y1 = np.zeros(Ndata1)\n",
    "X2 = np.random.multivariate_normal(mu2,cov2,size=Ndata2)\n",
    "Y2 = np.ones(Ndata2)\n",
    "\n",
    "X = np.concatenate((X1,X2))\n",
    "Y = np.concatenate((Y1,Y2)) \n",
    "\n",
    "\n",
    "data = [X,Y] # data set in the format used in the logistic regression model below\n",
    "\n",
    "\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "You can plot the data set using the code the following code"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "ndata = Y.shape[0]\n",
    "color_dict= {0:'red', 1 :'blue'}\n",
    "colors = [color_dict[Y[i]] for i in range(ndata)]\n",
    "plt.scatter(X[:,0],X[:,1], c=colors)\n",
    "plt.title('Data') \n",
    "plt.xlim([mu1[0]-3*np.sqrt(cov1[0,0]),+mu2[0]+3*np.sqrt(cov2[0,0])])\n",
    "plt.ylim([mu1[1]-3*np.sqrt(cov1[1,1]),+mu2[1]+3*np.sqrt(cov2[1,1])])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A suitable logistic regression model specifying a posterior distribution $p(\\beta | x,y)$ of a form as described above is implemented in the class \"BayesianLogisticRegression\". The implemented regression model assumes a Gaussian prior. The code below initializes the classifier model we will use in the following. \n",
    " - Note: In the code $q$ should be interpreted as $\\beta$, i.e., the output of our sampler corresponds to trajectory of samples distributed according to $p(\\beta | x,y)$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "model = models.BayesianLogisticRegression(prior_mean = np.zeros(data_dim),\n",
    "                                               prior_cov = 100*np.eye(data_dim),\n",
    "                                               data = [X,Y], \n",
    "                                         )"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "- Set up an Ensbemble Quasi Newton sampler to sample from the posterior of the Bayesion logistic regression model specified above. It is recommended to use 10 replicas and a stepsize of $\\Delta t = .01$. Ensure that the number of time steps is sufficiently large that you get good exploration of your sampling space. (Hint: You may want to inspect the trace plot of $\\beta_1$ and $\\beta_2$ (i.e., `op.traj_q[:,0]` and `op.traj_[:,1]`) for this)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "After you have obtained a sample trajectory using the Quasi Ensemble Newton method, you can ue the below code to plot the prediction $\\mathbb{P}[y_{pred} | x_{pred},x,y]$ for points on the grid $x_{pred} \\in$ grid. Where grid = xx1 $\\times$ xx2, with the vectors xx1 and xx2  as specified below. \n",
    "If computations take to much time you may want to adapt the values of Nx1, Nx2, Neval."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Nx1 = 50\n",
    "Nx2 = 50\n",
    "xx1 = np.linspace(mu1[0]-3*np.sqrt(cov1[0,0]),+mu2[0]+3*np.sqrt(cov2[0,0]),Nx1)\n",
    "xx2 = np.linspace(mu1[1]-3*np.sqrt(cov1[1,1]),+mu2[1]+3*np.sqrt(cov2[1,1]),Nx2)\n",
    "\n",
    "figb, axb = model.plot_prediction(op.traj_q, grid=[xx1,xx2], Neval=100 )\n",
    "plt.show()\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "If the above code does not work you may also use the below code fragment "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": false
   },
   "outputs": [],
   "source": [
    "Nx1 = 50\n",
    "Nx2 = 50\n",
    "xx1 = np.linspace(mu1[0]-3*np.sqrt(cov1[0,0]),+mu2[0]+3*np.sqrt(cov2[0,0]),Nx1)\n",
    "xx2 = np.linspace(mu1[1]-3*np.sqrt(cov1[1,1]),+mu2[1]+3*np.sqrt(cov2[1,1]),Nx2)\n",
    "Neval=100\n",
    "show_training_data = True \n",
    "\n",
    "q_trajectory = op.traj_q\n",
    "z= np.zeros([len(xx1),len(xx2)])\n",
    "modthin = q_trajectory.shape[0]//Neval\n",
    "for i in range(len(xx1)):\n",
    "    for j in  range(len(xx2)):\n",
    "        x = np.array([xx1[i],xx2[j]])\n",
    "        z[i,j] = model.predict_from_sample(q_trajectory[::modthin,:model.dim], x)[0]\n",
    "\n",
    "\n",
    "fig2, ax2 = plt.subplots()\n",
    "ax2.pcolor(xx1, xx2, z.transpose(), cmap='RdBu', vmin=0, vmax=1)\n",
    "cax = ax2.pcolor(xx1, xx2, z.transpose(), cmap='RdBu', vmin=0, vmax=1)\n",
    "cbar = fig2.colorbar(cax)\n",
    "cbar.ax.set_ylabel('$\\mathbb{P}(Y = 1)$')\n",
    "\n",
    "if show_training_data:\n",
    "    '''\n",
    "    Include training data\n",
    "    '''\n",
    "    X,Y = model.data\n",
    "    ndata = Y.shape[0]\n",
    "    color_dict= {0:'red', 1 :'blue'}\n",
    "    colors = [color_dict[Y[i]] for i in range(ndata)]\n",
    "\n",
    "    ax2.scatter(X[:,0],X[:,1], c=colors)\n",
    "\n",
    "ax2.set_title('Prediction')\n",
    "ax2.set_xlabel('$$')\n",
    "ax2.set_ylabel('$$')\n",
    "plt.show() "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "collapsed": true
   },
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "anaconda-cloud": {},
  "kernelspec": {
   "display_name": "Python [conda env:anaconda]",
   "language": "python",
   "name": "conda-env-anaconda-py"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.5.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}