09_Stochastic.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# 9. Stochastic Methods"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Some computations involve random numbers, such as simulating stochastic processes or searching in a high dimension space. \n",
    "\n",
    "References:\n",
    "* Christopher M. Bishop: Pattern Recognition and Machine Learning, Chapter 11: Sampling methods. Springer, 2006.\n",
    "* Jun S. Liu: Monte Carlo Strategies in Scientific Computing. Springer, 2004."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {
    "scrolled": true
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Random Numbers\n",
    "How can we generate *random* numbers from deterministic digital computers? They are, of course, *pseudo* random numbers."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The classic way of generating pseudo random numbers is *Linear Congruent Method* using a sequential dynimcs:  \n",
    "\n",
    "$$ x_{i+1}=ax_i+b \\mod m. $$\n",
    "\n",
    "Below is what was long used in the Unix system."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def lcm(n=1, seed=-1726662223):\n",
    "    \"\"\"Random integers <2**31 by linear congruent method\"\"\"\n",
    "    a = 1103515245\n",
    "    b = 12345\n",
    "    c = 0x7fffffff  # 2**31 - 1\n",
    "    x = np.zeros(n+1, dtype=int)\n",
    "    x[0] = seed\n",
    "    for i in range(n):\n",
    "        x[i+1] = (a*x[i] + b) & c  # bit-wise and\n",
    "    return(x[1:])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = lcm(1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1, 3, (1,2))\n",
    "plt.plot(x, '.')\n",
    "plt.ylabel('x')\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.hist(x, orientation='horizontal');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "A common problem with LCM is that lower digits can fall in simple cycles."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[:100]%4"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Also the samples can form crystal-like structure when consecutive numbers are taken as a vector."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "`numpy.random` uses an advanced method called *Mersenne Twister* which overcomes these problems."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.random.randint(2**31, size=1000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.subplot(1, 3, (1,2))\n",
    "plt.plot(x, '.')\n",
    "plt.ylabel('x')\n",
    "plt.subplot(1, 3, 3)\n",
    "plt.hist(x, orientation='horizontal');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x[:100]%4"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Integrating area/volume\n",
    "Uniform random numbers can be used to approximate integration\n",
    "\n",
    "$$ \\int_V f(x)dx \\simeq \\frac{|V|}{n}\\sum_{i=1}^n f(x_i) $$\n",
    "\n",
    "by uniform samples $x_i$ in a volume $V$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "For example, let us evaluate a volume of a sphere."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Sphere in 3D\n",
    "def sphere(x):\n",
    "    \"\"\"height of a half sphere\"\"\"\n",
    "    h2 = 1 - x[0]**2 - x[1]**2\n",
    "    return np.sqrt((h2>0)*h2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 10000\n",
    "x = np.random.random((2, m)) # one quadrant\n",
    "v = np.sum(sphere(x))/m\n",
    "print(8*v)\n",
    "print(4/3*np.pi)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Let us evaluate the volume of n-dimension sphere, which is supposed to be\n",
    "\n",
    "$$ \\frac{\\pi^{n/2}}{\\Gamma(n/2-1)} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "from scipy.special import gamma"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# n-dimensional sphere\n",
    "def nsphere(x):\n",
    "    \"\"\"height of a half sphere in n-dim\"\"\"\n",
    "    h2 = 1 - np.sum(x**2, axis=0)\n",
    "    return np.sqrt((h2>0)*h2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 10000\n",
    "for n in range(2,20):\n",
    "    x = np.random.random((n-1, m)) # one quadrant\n",
    "    v = np.sum(nsphere(x))/m\n",
    "    print(n, 2**n*v, np.pi**(n/2)/gamma(n/2+1))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Non-uniform Distributions\n",
    "How can we generate samples following a non-uniform distribution $p(x)$?\n",
    "\n",
    "If the cumulative density function \n",
    "\n",
    "$$f(x) = \\int_{-\\infty}^{x} p(u) du$$\n",
    "\n",
    "is known, we can map uniformly distributed samples $y_i\\in[0,1]$ to $x_i = f^{-1}(y_i)$."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Exponential distribution\n",
    "\n",
    "$$ p(x;\\mu) = \\frac{1}{\\mu}e^{-\\frac{x}{\\mu}} $$\n",
    "\n",
    "$$ f(x;\\mu) = 1 - e^{-\\frac{x}{\\mu}} $$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Exponential distribution in [0,infinity)\n",
    "def p_exp(x, mu=1):\n",
    "    \"\"\"density function of exponential distribution\"\"\"\n",
    "    return np.exp(-x/mu)/mu\n",
    "\n",
    "def f_exp(x, mu=1):\n",
    "    \"\"\"cumulative density function of exponential distribution\"\"\"\n",
    "    return 1 - np.exp(-x/mu)\n",
    "\n",
    "def finv_exp(y, mu=1):\n",
    "    \"\"\"inverse of cumulative density function of exponential distribution\"\"\"\n",
    "    return -mu*np.log(1 - y)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 10, 100)\n",
    "plt.plot(x, p_exp(x, 2))\n",
    "plt.plot(x, f_exp(x, 2))\n",
    "y = np.arange(0, 1, 0.01)\n",
    "plt.plot(finv_exp(y, 2), y, ':')\n",
    "plt.xlabel('x'); plt.ylabel('p, f');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def x_exp(n=1, mu=1):\n",
    "    \"\"\"sample from exponential distribution\"\"\"\n",
    "    ys = np.random.random(n)  # uniform in [0,1]\n",
    "    return finv_exp(ys, mu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x_exp(10, 2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = x_exp(1000, 4)\n",
    "plt.hist(xs, bins=20, density=True)\n",
    "plt.plot(x, p_exp(x, 4))\n",
    "plt.xlabel('x'); plt.ylabel('p');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "When a stochastic variable $x$ following the distribution $p(x)$ is transformed by $y=g(x)$, the distribution of $y$ is given by \n",
    "\n",
    "$$p(y) = \\left|\\frac{\\partial x}{\\partial y}\\right|p(x)$$\n",
    "\n",
    "where $|\\,|$ means the absolute value for a scalar derivative and the determinant for a Jacobian matrix."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Normal distribution\n",
    "A common way of generating a normal (Gaussian) distribution\n",
    "\n",
    "$$ p(x) = \\frac{1}{\\sqrt{2\\pi}}e^{-\\frac{1}{2}x^2} $$\n",
    "\n",
    "is to consider a 2D normal distribution\n",
    "\n",
    "$$ p(x_1,x_2) = p(x_1)p(x_2) = \\frac{1}{2\\pi}e^{-\\frac{1}{2}(x_1^2+x_2^2)}. $$\n",
    "\n",
    "For a polar coordinate tranformation\n",
    "\n",
    "$$ x_1 = r\\cos\\theta $$\n",
    "\n",
    "$$ x_2 = r\\sin\\theta, $$\n",
    "\n",
    "the Jacobian is \n",
    "\n",
    "$$ \\frac{\\partial (x_1,x_2)}{\\partial (r,\\theta)} = \n",
    "        \\left(\\begin{array}{cc} \\cos\\theta & -r\\sin\\theta\\\\\n",
    "        \\sin\\theta & r\\cos\\theta\\end{array}\\right)$$\n",
    "\n",
    "and its determinant is\n",
    "\n",
    "$$ \\det\\frac{\\partial (x_1,x_2)}{\\partial (r,\\theta)} = \n",
    "        r\\cos^2\\theta + r\\sin^2\\theta = r. $$\n",
    "\n",
    "Thus we have the relationship\n",
    "\n",
    "$$ p(r,\\theta) = r p(x_1,x_2) = \\frac{r}{2\\pi} e^{−\\frac{r^2}{2}}. $$\n",
    "\n",
    "By further transforming $u=r^2$, from $\\frac{du}{dr}= 2r$, we have \n",
    "\n",
    "$$ p(u,\\theta) = \\frac{1}{2r}p(r,\\theta) = \\frac{1}{4\\pi} e^{−\\frac{u}{2}}.$$\n",
    "\n",
    "Thus we can sample $u$ by exponential distribution \n",
    "\n",
    "$$p(u)=\\frac{1}{2} e^{−\\frac{u}{2}}$$\n",
    "\n",
    "and $\\theta$ by uniform distribution in $[0,2\\pi)$, and then transform them to $x_1$ and $x_2$ to generate two samples following normal distribution. \n",
    "\n",
    "This is known as Box-Muller method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def box_muller(n=1):\n",
    "    \"\"\"Generate 2n gaussian samples\"\"\"\n",
    "    u = x_exp(n, 2)\n",
    "    r = np.sqrt(u)\n",
    "    theta = 2*np.pi*np.random.random(n)\n",
    "    x1 = r*np.cos(theta)\n",
    "    x2 = r*np.sin(theta)\n",
    "    return np.hstack((x1,x2))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "box_muller(1)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "xs = box_muller(1000)\n",
    "plt.hist(xs, bins=20, density=True)\n",
    "# check how the histogram fits the pdf\n",
    "x = np.linspace(-5, 5, 100)\n",
    "px = np.exp(-x**2/2)/np.sqrt(2*np.pi)\n",
    "plt.plot(x, px)\n",
    "plt.xlabel('x'); plt.ylabel('p');"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In practice, you can sample by the whole variety of distributions in numpy.random package, but it's good to know the principle."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "#help(np.random)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Rejection Sampling\n",
    "What can we do when transformations from uniform distribution is not available?  \n",
    "\n",
    "Let us find a *proposal distribution* $q(x)$ for which samples are easily generated and covers the target distribution as $p(x)\\le cq(x)$ with a scaling constant $c>0$.  \n",
    "Then take a sample from $q(x)$ and accept it with probability $\\frac{p(x)}{cq(x)}$.\n",
    "\n",
    "### Gamma distribution\n",
    "*Gamma distribution* is an extension of exponential distribution defined as\n",
    "\n",
    "$$ p(x; k, \\theta) = \\frac{1}{\\Gamma(k)\\theta^k}x^{k-1}e^{-\\frac{x}{\\theta}} $$\n",
    "\n",
    "with the shape parameter $k>0$ and the scaling parameter $\\theta$. \n",
    "\n",
    "$\\Gamma(k)$ is the *gamma function*, which is a generalization of factorial, $\\Gamma(k)=(k-1)!$ for an integer $k$. \n",
    "\n",
    "Let us generate samples from gamma distribution with integer $k$ and $\\theta=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def p_gamma(x, k=1):\n",
    "    \"\"\"gamma distribution with integer k\"\"\"\n",
    "    return x**(k-1)*np.exp(-x)/np.prod(range(1,k))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = np.linspace(0, 10, 100)\n",
    "for k in range(1, 6):\n",
    "    plt.plot(x, p_gamma(x, k), label='k={}'.format(k))\n",
    "plt.xlabel('x'); plt.ylabel('p(x)'); plt.legend();"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Consider the exponential distribution $q(x;\\mu)=\\frac{1}{\\mu}e^{-\\frac{x}{\\mu}}$ as the proposal distribution.  \n",
    "For\n",
    "\n",
    "$$ \\frac{p(x;k)}{q(x;\\mu)} = \\frac{\\mu}{(k-1)!}x^{k-1}e^{-(1-\\frac{1}{\\mu})x}, $$\n",
    "\n",
    "we set\n",
    "\n",
    "$$ \\frac{d}{dx}\\frac{p(x;k)}{q(x;\\mu)} = 0 $$\n",
    "\n",
    "and have\n",
    "\n",
    "$$ \\left((k-1)x^{k-2}+\\frac{1-\\mu}{\\mu}x^{k-1}\\right)e^{\\frac{1-\\mu}{\\mu}x} = 0, $$\n",
    "\n",
    "$$ x = \\frac{\\mu(k-1)}{\\mu-1}, $$\n",
    "\n",
    "where $\\frac{p(x;k)}{q(x;\\mu)}$ takes the maximum\n",
    "\n",
    "$$ \\frac{\\mu^k}{(k-1)!}\\left(\\frac{k-1}{\\mu-1}\\right)^{k-1}e^{1-k}. $$\n",
    "\n",
    "By taking $\\mu=k$, \n",
    "\n",
    "$$ c=\\frac{k^k}{(k-1)!}e^{1-k} $$\n",
    "\n",
    "satisfies $p(x)\\le cq(x)$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "c = (k**k)/np.prod(range(1,k))*np.exp(1-k)\n",
    "#print(c)\n",
    "x = np.linspace(0, 10, 100)\n",
    "plt.plot(x, p_gamma(x, k))\n",
    "plt.plot(x, c*p_exp(x, k))\n",
    "plt.xlabel('x'); plt.ylabel('p, q');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def x_gamma(n=10, k=1):\n",
    "    \"\"\"sample from gamma distribution by rejection sampling\"\"\"\n",
    "    c = (k**k)/np.prod(range(1,k))*np.exp(1-k)\n",
    "    xe = x_exp(n, k)\n",
    "    paccept = p_gamma(xe, k)/(c*p_exp(xe, k))\n",
    "    accept = np.random.random(n)<paccept\n",
    "    xg = xe[accept]  # rejection sampling\n",
    "    #print(\"c =\", c, \"; acceptance rate =\", len(xg)/n)\n",
    "    return xg"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "k = 3\n",
    "xs = x_gamma(1000, k)\n",
    "print('accept:', len(xs))\n",
    "plt.hist(xs, bins=20, density=True)\n",
    "plt.plot(x, p_gamma(x, k))\n",
    "plt.xlabel('x'); plt.ylabel('p');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Importance Sampling\n",
    "One probelm of rejection sampling is that we may end up rejecting many of the samples when we cannot find a good proposal distribution.\n",
    "\n",
    "When the aim is not just to take samples from $p(x)$ but to take \n",
    "the *expectation* of a function $h(x)$ with respect to the distribution $p(x)$\n",
    "\n",
    "$$ E_p[h(x)] = \\int h(x)p(x)dx, $$\n",
    "\n",
    "we can use the ratio between the target and proposal distributions to better utilize the samples.\n",
    "\n",
    "$$ E_p[h(x)] = \\int h(x)p(x)dx = \\int h(x)\\frac{p(x)}{q(x)}q(x)dx\n",
    "    = E_q\\left[\\frac{p(x)}{q(x)}h(x)\\right] $$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "collapsed": true
   },
   "source": [
    "### Mean and variance\n",
    "Let us consider taking the mean\n",
    "\n",
    "$$ \\mu = E_p[x] = \\int xp(x)dx $$\n",
    "\n",
    "and the variance\n",
    "\n",
    "$$ \\sigma^2 = E_p[(x-\\mu)^2] = \\int (x-\\mu)^2p(x)dx  $$\n",
    "\n",
    "for the gamma distribution, which are known to be $k\\theta$ and $k\\theta^2$, respectively."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mv_exp(n=1, mu=1):\n",
    "    \"\"\"mean and variance of exponential distribution\"\"\"\n",
    "    x = x_exp(n, mu)\n",
    "    mean = np.mean(x)\n",
    "    var = np.var(x)\n",
    "    return (mean, var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv_exp(1000, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mv_gamma(n=1, k=1):\n",
    "    \"\"\"mean and variance of gamma distribution by rejection sampling\"\"\"\n",
    "    x = x_gamma(n, k)  # by rejection sampling\n",
    "    mean = np.mean(x)\n",
    "    var = np.var(x)\n",
    "    return (mean, var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv_gamma(100, 3)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def mv_gamma_is(n=1, k=1):\n",
    "    \"\"\"mean and variance of gamma distribution by importance sampling\"\"\"\n",
    "    x = x_exp(n, k)  # sample by exponential distribution\n",
    "    importance = p_gamma(x, k)/p_exp(x, k)\n",
    "    mean = np.dot(importance, x)/n\n",
    "    var = np.dot(importance, (x-mean)**2)/(n - 1)\n",
    "    return (mean, var)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "mv_gamma_is(1000, 3)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Compare the variability of the estimate of the mean."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "m = 100  # number of runs\n",
    "n = 100  # number of samples\n",
    "k = 3\n",
    "means = np.zeros((m, 2))\n",
    "for i in range(m):\n",
    "    means[i,0], var = mv_gamma(n, k)\n",
    "    means[i,1], var = mv_gamma_is(n, k)\n",
    "print(\"RS: \", np.mean(means[:,0]), np.var(means[:,0]))\n",
    "print(\"IS: \", np.mean(means[:,1]), np.var(means[:,1]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Markov Chain Monte Carlo (MCMC)\n",
    "What if we don't know the right proposal distribution? We can take a nearby point of a previous sample as a candidate sample and keep it if the candidate's probability is high. The method is callel Markov Chain Monte Carlo (MCMC) method and used practically for sampling from unknow distributions in high dimensional space."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Metropolis Sampling\n",
    "A simple example of MCMC is Metropolis sampling, which requires only unnormalized propability $\\tilde{p}(x)\\propto p(x)$ of samples for relative comparison.\n",
    "\n",
    "A new candidate $x^*$ is generated by a symmetric proposal distribution $q(x^*|x_k)=q(x_k|x^*)$, such as a gaussian distribution, and acctepted with the probability of\n",
    "\n",
    "$$ p_\\mbox{acc} = \\min(1, \\frac{p(x^*)}{p(x_k)}) $$ "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def metropolis(p, x0, sig=0.1, m=1000):\n",
    "    \"\"\"metropolis: Metropolis sampling\n",
    "    p:unnormalized probability, x0:initial point,\n",
    "    sif:sd of proposal distribution, m:number of sampling\"\"\"\n",
    "    n = len(x0)  # dimension\n",
    "    p0 = p(x0)\n",
    "    x = []\n",
    "    for i in range(m):\n",
    "        x1 = x0 + sig*np.random.randn(n)\n",
    "        p1 = p(x1)\n",
    "        pacc = min(1, p1/p0)\n",
    "        if np.random.rand()<pacc:\n",
    "            x.append(x1)\n",
    "            x0 = x1\n",
    "            p0 = p1\n",
    "    return np.array(x)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def croissant(x):\n",
    "    \"\"\"croissant-like distribution in 2D\"\"\"\n",
    "    return np.exp(-x[0]**2 - (x[1]-x[0]**2)**2)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "r = 4  # plot rage\n",
    "x = np.linspace(-r, r)\n",
    "X, Y = np.meshgrid(x, x)\n",
    "Z = croissant(np.array([X,Y]))\n",
    "plt.contour(X, Y, Z);"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = metropolis(croissant, [0,0])\n",
    "s = len(x); print(s)\n",
    "plt.contour(X, Y, Z)\n",
    "plt.scatter(x[:,0], x[:,1], c=np.arange(s), marker='.');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "x = metropolis(croissant, [2,0], sig=0.1, m=1000)\n",
    "s = len(x); print(s)\n",
    "plt.contour(X, Y, Z)\n",
    "plt.scatter(x[:,0], x[:,1], c=np.arange(s), marker='.');"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "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.12.2"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}