Differences

This shows you the differences between two versions of the page.

--- cs401r_w2016:lab10 [2016/03/17 19:14]
admin
+++ cs401r_w2016:lab10 [2021/06/30 23:42] (current)
@@ Line 23: / Line 23: @@
 {{ :cs401r_w2016:lab9_t2.png?direct&400 |}}
 {{ :cs401r_w2016:lab9_t3.png?direct&400 |}}
+Your notebook should produce similar plots for the HMC algorithm, although you only need to produce two plots (one histogram, and one state evolution plot, instead of three of each). So to clarify, your notebook should have 6 plots for part one (three histograms, three evolution plots) and two plots for part two (one histogram, one evolution plot).
+Your notebook should also include a small writeup of your results, as described below.
 ----
@@ Line 30: / Line 34: @@
   * 25% Correct implementation of Metropolis Hastings inference
-  * 25% Correct calculation of gradients
+  * 5% Correct calculation of gradients
-  * 25% Correct implementation of Hamiltonian MCMC
+  * 45% Correct implementation of Hamiltonian MCMC
-  * 25% Final plot is tidy and legible
+  * 15% A small write-up comparing and contrasting MH, HMC, and the different proposal distributions
+  * 10% Final plot is tidy and legible
 ----
 ====Description:====
-For this lab, you will implement the basic particle filter found in Algorithm 23.1 of the book (pg. 828).  You can assume that your particle proposal distribution is equal to the dynamics prior (see Section 23.5.4), which leads to a very simple form of particle filtering.  Specifically, you may assume that the "dynamics" are that the moves in a random direction and changes angles randomly.  I used $p(z_{t+1}|z_{t}) = \mathcal{N}(0,\Sigma)$, where $\Sigma$ is a parameter you will need to tune.
+For this lab, you will code two different MCMC algorithms.  Each will attempt to draw samples from the same distribution, given by the following density function:
-You will need the following files:
+<code python>
+import numpy as np
+def p( x, temperature=1.0 ):
+    return np.exp( -10*temperature*((x-2)**2) ) + 0.3*np.exp( -0.5*10*temperature*((x+1)**2) )
+</code>
-[[http://hatch.cs.byu.edu/courses/stat_ml/sensors.mat|sensors.mat]] contains the data you will use (you can also regenerate it directly, using the script below).  The data is in the ''sonars'' field, and the true states are in the ''true_states'' field.  ''true_states'' is a 3xT matrix, where each 3x1 column represents the true x,y and angle (\in 0,2*pi) of the robot at time t.
+This distribution has two modes that are separated by a region of low probability.
-[[http://hatch.cs.byu.edu/courses/stat_ml/lab8_common.py|lab8_common.py]] contains utility functions that will enable you to create a simple map, and compute what a robot would see from a given position and angle.
+**Part 1: Metropolis Hastings**
-[[http://hatch.cs.byu.edu/courses/stat_ml/lab8_gen_sonars.py|lab8_gen_sonars.py]] is the script that I used to generate the data for this lab.  If you run it, it will generate a new dataset, but it will also produce a nice visualization of the robot as it moves around, along with what the sensors measure.
+For this part, you should implement the basic MH MCMC algorithm.  You should use a Gaussian proposal distribution with three different variances ''(0.1, 1.0 and 10.0)''.  Your sampler should start with an initial state of ''0''.
-To get started, check out the provided functions in lab8_common.py.  There are two critical functions:
+For each different proposal distribution, you should run your MCMC chain for 10,000 steps, and record the sequence of states.  Then, you should produce a visualization of the distribution of states, and overlay a plot of the actual target distribution.  They may or may not match (see, for example, the first example plot in the Description section).
-  - ''create_map'' will create a map object and return it
+Furthermore, for each proposal distribution, you should run three independent chains (you can do these sequentially or in parallel, as you like).  You should display each of these three chains on a single plot with time on the x-axis and the state on the y-axis.  Ideally, you will see each of the three chains mixing between two modes; you may notice other features of the behavior of the samplers as well, which you should report in your writeup!
-  - ''cast_rays'' will simulate what a robot would see from a specific x,y,theta position, given a map
-''cast_rays'' is fully vectorized, meaning that you can simulate from multiple locations simultaneously.  For example, in my code, I only call ''cast_rays'' once:
+**Part 2: Hamiltonian MCMC**
+For this part, you will code the Hamiltonian MCMC algorithm, as discussed in class. You will run three independent chains and report them in the same graphs.  To do this, you will need to compute the gradient of the density function with respect to the state.  An easy easy way to do this is to use the [[https://github.com/HIPS/autograd|autograd]] package:
 <code python>
-    part_sensor_data = cast_rays( parts, room_map )
+from autograd import grad
-    liks = ray_lik( data[:,t:t+1], part_sensor_data.T ).T
+import autograd.numpy as np
+grad_p = grad( p )
 </code>
+The function ''grad_p'' accepts the same parameters as ''p'', but it returns the gradient, instead of the density.
-where I have represented my particles as a 3xN matrix.  In the above example call, ''part_sensor_data'' is an Nx11 matrix.  ''ray_lik'' is my implementation of the likelihood function $p(y_t|z_t)$, and is also fully vectorized.  It returns an Nx1 matrix, with the likelihood of each particle.
+You should use the leapfrog method to integrate the dynamics.
-Your proposal distribution can be anything you want.  Be creative.  You may have to experiment!
+Remember that you will need to introduce as many momentum variables as there are state (ie, position) variables.
-You should experiment with two different likelihoods.  First, try using a simple Gaussian likelihood.  You will need to tune the variance, but don't spend too long on it; I never got this to work very well.  Second, try using a simple Laplacian distribution.  (The Laplacian is like a Gaussian, except that instead of squaring the difference, one uses an absolute value): $$p(y_t|z_t^s) = \exp \left( -\sigma \sum_i \lvert (y_t)_i - \mathrm{cast\_rays}(z_t^s)_i \rvert \right) $$, where $\sigma$ is a parameter that you will need to tune.
+A detailed explanation of Hamiltonian MCMC can be found here:[[http://www.mcmchandbook.net/HandbookChapter5.pdf|Hamiltonian MCMC]].
-You will have to experiment with the parameters of these distributions to make them work.  Try values between 0.01 and 1.
+  * You will find the equations describing the leapfrog method in Equations 5.18, 5.19 and 5.20.
+  * You will find a description of how to convert a given ''p(x)'' into a Hamiltonian in Section 5.3.1.
+  * You will find a description of the complete HMC algorithm in section 5.3.2.1
-----
+Remember that you will alternate between two steps:
-====Hints:====
-You may find the ''numpy.random.mtrand.multinomial'', or ''numpy.random.choice'' functions helpful.
+  - Sampling the momentum conditioned on the position.  This is just sampling from a Gaussian.
+  - Proposing a new state for the position, given the momentum.  This involves integrating the dynamics, and then accepting or rejecting based on integration error.
-It can be frustrating to experiment with a particle filter when the interface is too slow.  By using ''pyqtgraph'' instead of ''matplotlib'', you can speed things up considerably.
+You will have to tune two parameters in order to implement HMC: the variance of the momentum variables, and the timestep used for integrating the dynamics.  Experiment with both, and report your results using plots like those you prepared for Part 1.
-Here is some code to create a ''pyqtgraph'' window, and show the map.
-<code python>
+**Part 3: Observations**
-from pyqtgraph.Qt import QtGui, QtCore
-import pyqtgraph as pg
-app = QtGui.QApplication( [] )
+You have now coded two different inference algorithms, and a few variants of each.  For this section, you must provide a small write-up that compares and contrasts each.  Answer at least the following questions:
-win = pg.GraphicsWindow( title="Particle filter" )
-win.resize( 600, 600 )
-win.setWindowTitle( 'Particle filter' )
-pg.setConfigOptions( antialias=True )
-p3 = win.addPlot( title="Room map" )
+  - What was the acceptance rate of each algorithm?  (ie, what percentage of proposals were accepted)
+  - Why don't some inference algorithms explore both modes of the density?
+  - Why do some algorithms stay in the same state repeatedly?  Is this good or bad?
+  - What were the best values for the variance of the momentum variables and the timestep you found?  How did you know that they were good?
-for i in range( 0, room_map.shape[0] ):
+----
-    p3.plot( [room_map[i,0], room_map[i,2]], [room_map[i,1], room_map[i,3]] )
+====Hints:====
-p3.setXRange( -0.1, 1.1 )
-p3.setYRange( -0.1, 1.1 )
-</code>
-To speed things up further, you should use the ''.setData'' method on your plots.  This replaces the data used for a particular plotting element, and is faster than deleting and recreating the plot element.  This is particularly important when trying to update many particles.
+You may find ''plt.hist'' with the ''normed=True'' option helpful.
-Here's a snippet of my code.  ''ts_plot'' and ''ex_plot'' are initialized to ''None''.  My particles are stored in a 3xN matrix called ''parts'', my expected states are stored in a 3xt matrix called ''exs''.
-<code python>
-if ts_plot == None:
-    ts_plot = p3.plot( true_states[0,0:t+1], true_states[1,0:t+1], pen=(0,0,255) )
-    ex_plot = p3.plot( exs[0,0:t+1], exs[1,0:t+1], pen=(0,255,0) )
-    pts = p3.scatterPlot( parts[0,:], parts[1,:], symbol='o', size=1, pen=(255,100,100) )
-else:
-    ts_plot.setData( true_states[0,0:t+1], true_states[1,0:t+1] )
-    ex_plot.setData( exs[0,0:t+1], exs[1,0:t+1] )
-    pts.setData( parts[0,:], parts[1,:] )
-pg.QtGui.QApplication.processEvents()
-</code>
-(you must call ''processEvents'' to actually update the display)
-The combination of these two things allowed me to experiment with my particle filter in real-time.  Note that you are not required to produce any sort of visualization over time, just a single plot at the end of running your code.  The use of ''pyqtgraph'' is strictly for the work you will do debugging your particle filter.

BYU CS classes

User Tools

Site Tools

Differences

Page Tools