{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# Sampling"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"###### COMP4670/8600 - Statistical Machine Learning - Tutorial\n",
"\n",
"## Textbook Questions\n",
"These questions are hand picked to both be of reasonable difficulty and demonstrate what you are expected to be able to solve. The questions are labelled in Bishop as either $\\star$, $\\star\\star$, or $\\star\\star\\star$ to rate its difficulty.\n",
"\n",
"- **8.1**\n",
"- **8.2**\n",
"- **8.10** Hint: For the first part try marginalise the joint distribution to work out p(a,b). For the second part you may want to use Bayes' rule.\n",
"- **8.13**\n",
"- **8.14** Hint: How does energy relate to the distribution of x and y?\n",
"- **8.15** Hint: You can follow the same steps used in the single variable example at the start of 8.4.1\n",
"- **8.18** (Challenge) Hint: To convert to an undirected graph look at the start of 8.3.4. To convert to a directed graph simply pick a node to be your root node.\n",
"\n",
"### Assumed knowledge\n",
"- Sampling (lectures)\n",
"\n",
"### After this lab, you should be comfortable with:\n",
"- The concept of Monte Carlo methods in general\n",
"- Importance sampling in particular"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Setting up the environment\n",
"$\\newcommand{\\Ex}{\\mathbb{E}}$\n",
"$\\newcommand{\\dd}{\\mathrm{d}}$\n",
"$\\newcommand{\\DUniform}[3]{\\mathscr{U}\\left(#1 ~\\middle|~ #2, #3\\right)}$\n",
"$\\newcommand{\\DNorm}[3]{\\mathscr{N}\\left(#1 ~\\middle|~ #2, #3\\right)}$\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:28:39.790440Z",
"start_time": "2019-05-12T12:28:39.526916Z"
}
},
"outputs": [],
"source": [
"import math\n",
"import numpy as np\n",
"from scipy.stats import uniform, multivariate_normal\n",
"import matplotlib.pyplot as plt\n",
"\n",
"%matplotlib inline"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"The aim of this tutorial is to investigate the effects of different proposal distributions on importance sampling.\n",
"\n",
"## Sampling from Gaussian and Uniform Distributions\n",
"\n",
"*Note that help for a certain function can be obtained by using the question mark, for example * ```?uniform``` or ```?norm```.\n",
"\n",
"Repeat the following twice, once each for Gaussian with zero mean and unit variance and Uniform on the unit square.\n",
"1. Sample 1000 data points from a two dimensional distribution.\n",
"2. Compute the [two dimensional histogram](http://docs.scipy.org/doc/numpy/reference/generated/numpy.histogram2d.html), with 10 bins in each dimension.\n",
"3. Visualise the histogram as a heatmap. There are [various ways of doing this](http://thomas-cokelaer.info/blog/2014/05/matplotlib-difference-between-pcolor-pcolormesh-and-imshow/) as a two dimensional image. Use an [appropriate colormap](https://jakevdp.github.io/blog/2014/10/16/how-bad-is-your-colormap/).\n",
"4. Visualise the difference between the theoretical and empirical values of the density.\n",
"\n",
"The aim of this exercise is to observe the challenges of sampling in more than 1 dimension."
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:28:40.063429Z",
"start_time": "2019-05-12T12:28:39.791843Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:30:10.854362Z",
"start_time": "2019-05-12T12:30:10.571147Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## A function for performing estimation\n",
"\n",
"Given are the following function $f(x, y)$\n",
"\n",
"\\begin{equation*}\n",
" f(x, y) = x \\, y \\cos(x) \\, \\cos(y),\n",
"\\end{equation*}\n",
"\n",
"and the (unnormalised!) distribution $ \\widetilde{p}(x, y) $\n",
"\\begin{equation*}\n",
" \\widetilde{p}(x, y) = \n",
" \\exp \\left\\{- \\frac{1}{4}((x - 2)^2 + (y - 3)^2) \\right\\}\n",
" - \\frac{3}{4} \\exp \\left\\{- \\frac{1}{2} ((x - 2)^2 + (y - 3)^2) \\right\\}\n",
"\\end{equation*}\n",
"\n",
"The goal is to numerically estimate $ \\Ex_{p(x, y)}[f(x, y)] $ defined as\n",
"\\begin{equation*}\n",
" \\Ex_{p(x, y)}[f(x, y)] = \\int_{-\\infty}^{\\infty} f(x, y) \\, p(x, y) \\dd x \\dd y\n",
"\\end{equation*}\n",
"\n",
"where $ p(x, y) $ is the normalised probability distribution proportional to $\\widetilde{p}(x, y) $, i.e.\n",
"\\begin{equation*}\n",
" p(x, y) = \\frac{ \\widetilde{p}(x, y) }\n",
" { \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} \\widetilde{p}(x, y) \\dd x \\dd y }\n",
"\\end{equation*}\n",
"\n",
"Note, that with importance sampling, this can be achieved without calculating the normalisation. All we need to do, is to use sample points from an appropriate distribution $ \\widetilde{q}(x, y) $ \n",
"which has most of its probability mass in regions where $ \\widetilde{p}(x, y) $ is also nonzero.\n",
"\n",
"In the following, we will implement importance sampling to estimate $ \\Ex_{p(x, y)}[f(x, y)] $ for two choices of $ \\widetilde{q}(x, y) $, a Uniform\n",
"Distribution and a Gaussian Distribution.\n",
"\n",
"**(optional) Plot the functions $f(x,y)$ and $\\widetilde{p}(x, y)$. Look at the 3D plotting section in [this tutorial](\\widetilde{p}(x, y)).**"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:28:40.754770Z",
"start_time": "2019-05-12T12:28:40.345956Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importance Sampling using the Uniform Distribution\n",
"\n",
"Draw random samples $ x_n $ and $ y_n $, $ n = 1, \\dots, N $, \n",
"each i.i.d. from the uniform distribution on the interval $[-10, 20] $. \n",
"\n",
"Now use the samples from the distribution $ \\widetilde{q}(x, y) $ to estimate \n",
"$ \\Ex_{p(x, y)}[f(x, y)] $ via importance sampling.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:28:59.816073Z",
"start_time": "2019-05-12T12:28:40.755936Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Importance Sampling using the Gaussian Distribution\n",
"\n",
"create random vectors from the following Gaussian Distribution\n",
"$ \\widetilde{q}(x, y) = \\DNorm{(x, y)^T}{\\mathbf{\\mu} = (2, 3)^T}{\\mathbf{\\Sigma} = \\begin{bmatrix} 2 & 0 \\\\ 0 & 3 \\\\ \\end{bmatrix}} $.\n",
"\n",
"Now use the samples from the distribution $ \\widetilde{q}(x, y) $ to estimate \n",
"$ \\Ex[p(x, y)]{f(x, y)} $ via importance sampling.\n",
"\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:29:12.341438Z",
"start_time": "2019-05-12T12:28:59.817235Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Comparing the empirical and analytic results\n",
"\n",
"For the function $f(x, y)$ and the normalised distribution $p(x, y)$,\n",
"the correct result can be computed analytically\n",
"\\begin{align*}\n",
" \\Ex_{p(x, y)}[f(x, y)]\n",
" & = \\int_{-\\infty}^{\\infty} \\int_{-\\infty}^{\\infty} f(x, y) \\, p(x, y) \n",
" \\dd x \\dd y \\\\\n",
" & = \\frac{16 (5 \\cos(1) + \\cos(5) + \\sin(1) - 5 \\sin(5))\n",
" - 3 e (7 \\cos(1) + 5 \\cos(5) + \\sin(1) - 5 \\sin(5)) }\n",
" {10 \\, e^2} \\\\ \n",
" & \\approx 0.670859\n",
"\\end{align*}\n",
"\n",
"Compare the convergence rate of the two approaches given above. Plot on the same plot the two curves showing the empirical expectation as a function of the number of samples as well as the analytical value.\n"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {
"ExecuteTime": {
"end_time": "2019-05-12T12:29:12.463468Z",
"start_time": "2019-05-12T12:29:12.342639Z"
}
},
"outputs": [],
"source": [
"# replace this with your solution, add and remove code and markdown cells as appropriate"
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"Discuss how one could choose an appropriate proposal distribution for a particular function $f(x,y)$ and distribution $p(x,y)$."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Answer\n",
"--- replace this with your solution, add and remove code and markdown cells as appropriate ---"
]
},
{
"cell_type": "markdown",
"metadata": {
"collapsed": true
},
"source": [
"## Interpretation\n",
" \n",
"1. How this setup can be used to estimate a posterior predictive mean (very easy), and a posterior predictive variance (slightly trickier).\n",
"- An example case in Bayesian machine learning which would would require sampling; i.e. be analytically intractable."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"### Answer\n",
"--- replace this with your solution, add and remove code and markdown cells as appropriate ---"
]
},
{
"cell_type": "code",
"execution_count": null,
"metadata": {},
"outputs": [],
"source": []
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3",
"language": "python",
"name": "python3"
},
"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.7.3"
}
},
"nbformat": 4,
"nbformat_minor": 1
}