accept-reject-method.ipynb

{
 "cells": [
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Simulation methods: laws using accept-reject method\n",
    "\n",
    "We here implement functions that generate samples following a specific law, from the random number generator.\n",
    "This implementation will use the accept-reject method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "import matplotlib.pyplot as plt\n",
    "from math import factorial\n",
    "from scipy.integrate import quad\n",
    "%matplotlib inline"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The `quad` function from the `scipy.integrate` module is used to compute the numerical integral of the law."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Absolute value normal law\n",
    "\n",
    "We first implement the absolute value normal law, which is defined as follows:\n",
    "$$f(x)=\\frac{2}{\\sqrt{2\\pi}}\\exp\\left(-\\frac{x^2}{2}\\right)\\ \\text{when}\\ x>0\\ \\text{else}\\ 0$$\n",
    "As we don't want to generate negative values, we will use the accept-reject method.\n",
    "We notice that:\n",
    "$$\\forall x>0,\\qquad f(x)\\leq\\sqrt{\\frac{2\\mathrm e}\\pi}\\mathrm{e}^{-x}$$\n",
    "Therefore, we set $c=\\sqrt{\\frac{2\\mathrm e}\\pi}$ and we generate $Y$ following an exponential law with $\\lambda=1$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def exponentielle(lam, n):\n",
    "    L=[]\n",
    "    for i in range(n):\n",
    "        u = np.random.random()\n",
    "        L.append(-np.log(1-u)/lam)\n",
    "    return L\n",
    "\n",
    "def abs_norm(n):\n",
    "    c = np.sqrt(2 * np.exp(1) / np.pi)\n",
    "    L=[]\n",
    "    while len(L) < n:\n",
    "        u = np.random.random()\n",
    "        g = exponentielle(1, 1)[0]\n",
    "        f= 2 / np.sqrt(2 * np.pi) * np.exp(-g**2/2)\n",
    "        h = f / (c * g)\n",
    "        if u <= h:\n",
    "            L.append(g)\n",
    "    return L\n",
    "\n",
    "X=abs_norm(10000)\n",
    "plt.hist(X, density= True, ec = 'white')\n",
    "f = lambda x: 2 / np.sqrt(2 * np.pi) * np.exp(-x**2/2)\n",
    "plt.plot(np.linspace(0, 3, 100), f(np.linspace(0, 3, 100)))\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Law on $\\mathbb R^2$\n",
    "\n",
    "We now implement a law on $\\mathbb R^2$, more specially around the unit circle.\n",
    "For that, we will use the following algorithm:\n",
    "1. Do `u1 <- 2 * rand - 1` and `u2 <- 2 * rand - 1` while `u1^2 + u2^2 > 1`.\n",
    "2. `U1 <- u1` and `U2 <- u2`.\n",
    "\n",
    "Then we plot the result, and also the unit circle, and a sqare of side 2 centered on the origin."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def va(n = 10000):\n",
    "    X=[]\n",
    "    Y=[]\n",
    "    for i in range(n):\n",
    "        u1 = 2 * np.random.random() - 1\n",
    "        u2 = 2 * np.random.random() - 1\n",
    "        while u1 * u1 + u2 * u2 <= 1:\n",
    "            u1 = 2 * np.random.random() - 1\n",
    "            u2 = 2 * np.random.random() - 1\n",
    "        X.append(u1)\n",
    "        Y.append(u2)\n",
    "    return X, Y\n",
    "\n",
    "sample1, sample2 = va(1000)\n",
    "u=0\n",
    "v=0\n",
    "a=1\n",
    "b=1\n",
    "t = np.linspace(0, 2*np.pi, 100)\n",
    "plt.figure(figsize=(10, 10))\n",
    "plt.scatter(sample1, sample2)\n",
    "plt.plot(u + a * np.cos(t), v + b * np.sin(t), color='red')\n",
    "plt.plot([-1, 1, 1 ,-1, -1], [1, 1, -1, -1, 1], color='green')\n",
    "plt.axis('equal')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Law on $\\mathbb R^3$\n",
    "\n",
    "We now consider a law on $\\mathbb R^3$, more specially around the unit ball.\n",
    "We will use a samilar algorithm as the previous one."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def va_boule(n):\n",
    "    L=[]\n",
    "    while len(L) < n:\n",
    "        x = -1 + 2 * np.random.random(3)\n",
    "        if np.linalg.norm(x) <= 1:\n",
    "            L.append(x)\n",
    "    return L\n",
    "\n",
    "sample = va_boule(1000)\n",
    "ax = plt.figure().add_subplot(projection='3d')\n",
    "ax.scatter([x[0] for x in sample], [x[1] for x in sample], [x[2] for x in sample])\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Gaussian binomial law\n",
    "\n",
    "### Introduction\n",
    "\n",
    "Using the accept-reject method, we can generate samples following a gaussian binomial law.\n",
    "To do so, we will use the Box-Muller method to generate samples following a gaussian law."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian2(n):\n",
    "    r = np.array(exponentielle(0.5, n))\n",
    "    u = 2 * np.pi * np.random.random(n)\n",
    "    return np.sqrt(r) * np.cos(u), np.sqrt(r) * np.sin(u)"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Law on $\\mathbb R^2$\n",
    "\n",
    "We now consider a normal law on $\\mathbb R^2$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "X,Y = gaussian2(10000)\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.scatter(X,Y)\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Law on $\\mathbb R$\n",
    "\n",
    "We are going to use this method to generate a normal sample of size 10000, and then we plot the histogram of the result.\n",
    "Thus, we can compare that with the normal law."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "gauss1 = gaussian2(10000)[0]\n",
    "plt.figure(figsize=(10, 7))\n",
    "plt.hist(gauss1, density= True, ec = 'white',label='Sample')\n",
    "plt.plot(np.linspace(min(gauss1), max(gauss1), 100), 1/np.sqrt(2*np.pi) * np.exp(-np.linspace(min(gauss1), max(gauss1), 100)**2/2),label='Gaussian law')\n",
    "plt.legend()\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We notice that the histogram is not exactly the same as the normal law, but it is close enough."
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Gaussian binomial law on $\\mathbb R^2$\n",
    "\n",
    "Now, we are going to simulate a Gaussian binomial law.\n",
    "For that, we set:\n",
    "1. $Z_1=\\mu_1+\\sigma_{1,1}X_1+\\sigma_{1,2}X_2$.\n",
    "2. $Z_2=\\mu_2+\\sigma_{2,1}X_1+\\sigma_{2,2}X_2$.\n",
    "\n",
    "where $X_1$ and $X_2$ are two independent normal samples.\n",
    "Then we define:\n",
    "$$\\mu=\\binom{\\mu_1}{\\mu_2}$$\n",
    "and\n",
    "$$\\sigma=\\begin{pmatrix}\\sigma_1^2&\\rho\\sigma_1\\sigma_2\\\\\\rho\\sigma_1\\sigma_2&\\sigma_2^2\\end{pmatrix}$$\n",
    "with:\n",
    "1. $\\sigma_1=\\sqrt{\\sigma_{1,1}^2+\\sigma_{1,2}^2}$.\n",
    "2. $\\sigma_2=\\sqrt{\\sigma_{2,1}^2+\\sigma_{2,2}^2}$.\n",
    "3. $\\rho=\\frac{\\sigma_{1,1}\\sigma_{2,1}+\\sigma_{1,2}\\sigma_{2,2}}{\\sigma_1\\sigma_2}$.\n",
    "\n",
    "The law defined by $Z=(Z_1,Z_2)$ is a Gaussian binomial law with parameters $\\mu$ and $\\sigma$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "def gaussian_binomial_density(z1, z2, sigma):\n",
    "    sigma1 = np.linalg.norm(sigma1)\n",
    "    sigma2 = np.linalg.norm(sigma2)\n",
    "    rho = (sigma[1,1] * sigma[1,2] + sigma[2,1] * sigma[2,2]) / (sigma1 * sigma2)\n",
    "    z = z1 ** 2 / sigma1 ** 2 + z2 ** 2 / sigma2 ** 2 - 2 * rho * z1 * z2 / (sigma1 * sigma2)\n",
    "    return 1 / (2 * np.pi * sigma1 * sigma2 * np.sqrt(1 - rho ** 2)) * np.exp(-z / (2 * (1 - rho ** 2)))"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In our case, we take $\\mu_1=0$ and $\\mu_2=0$.\n",
    "We also want to have $\\rho$ that takes the values $0.1$, $0.5$ and $0.9$.\n",
    "We need to find values for the $\\sigma_{i,j}$ in accordance with that.\n",
    "To do so, we set $\\sigma_{1,1}=1$ and $\\sigma_{2,2}=1$.\n",
    "Then we have:\n",
    "1. $\\sigma_{1,2}=0.01$.\n",
    "2. $\\sigma_{2,1}=0.09$.\n",
    "\n",
    "In this situation, we have $\\rho=0.1$.\n",
    "\n",
    "To have $\\rho = 0.5$, we set:\n",
    "1. $\\sigma_{1,2}=0.1$.\n",
    "2. $\\sigma_{2.1}=0.45$.\n",
    "\n",
    "To have $\\rho = 0.9$, we set:\n",
    "1. $\\sigma_{1,2}=0.5$.\n",
    "2. $\\sigma_{2.1}=0.77$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "sigma11 = 1\n",
    "sigma22 = 2\n",
    "sigma_list = [(0.01,0.1), (0.1,0.42), (0.5,0.55)]\n",
    "\n",
    "for s in sigma_list:\n",
    "    sigma12, sigma21 = s\n",
    "    sigma1 = np.linalg.norm([sigma11, sigma12])\n",
    "    sigma2 = np.linalg.norm([sigma21, sigma22])\n",
    "    rho = (sigma12 * sigma22 + sigma21 * sigma22) / (sigma1 * sigma2)\n",
    "    sigma = np.array([[sigma1 ** 2, rho * sigma1 * sigma2], [rho * sigma1 * sigma2, sigma2 ** 2]])\n",
    "    print(f'Computed rho = {rho}')\n",
    "    print(f'Computed sigma = {sigma}')\n",
    "    X1,X2 = gaussian2(10000)\n",
    "    Z1 = sigma11 * X1 + sigma12 * X2\n",
    "    Z2 = sigma21 * X1 + sigma22 * X2\n",
    "    plt.title(f'rho = {rho}')\n",
    "    plt.scatter(range(10000), Z2, label='Z2', color = 'blue')\n",
    "    plt.scatter(range(10000), Z1, label='Z1', color = 'red')\n",
    "    plt.legend()\n",
    "    plt.show()\n",
    "    plt.title(f'rho = {rho}, Z1')\n",
    "    plt.hist(Z1, density= True, ec = 'white')\n",
    "    plt.show()\n",
    "    plt.title(f'rho = {rho}, Z2')\n",
    "    plt.hist(Z2, density= True, ec = 'white')\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## Mixture of laws\n",
    "\n",
    "We are looking for a way to represent a mixture of law.\n",
    "Hence, we consider $f$ defined as follow:\n",
    "$$f(x)=p_1f_1(x)+p_2f_2(x)+p_3f_3(x)$$\n",
    "where $p_1=\\frac16$, $p_2=\\frac13$, $ p_3=\\frac12$, and $f_i$ are laws.\n",
    "In our case, we have:\n",
    "1. $f_1(x)=1$ if $x\\in[0,1]$ else $0$.\n",
    "2. $f_2(x)=\\frac12(2x-1)$ if $x\\in]1,2]$ else $0$.\n",
    "3. $f_3(x)=\\frac23(-3x+9)$ if $x\\in]2,3]$ else $0$.\n",
    "\n",
    "We can first plot the function $f$, between $0$ and $3$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "P = [1/6,1/3,1/2]\n",
    "\n",
    "f1 = lambda x : 1 if 0 <= x <= 1 else 0\n",
    "f2 = lambda x : 1/2*(2*x-1) if 1 < x <= 2 else 0\n",
    "f3 = lambda x : 2/3*(-3*x+9) if 2 < x <= 3 else 0\n",
    "\n",
    "def f(x):\n",
    "    A=np.zeros(len(x))\n",
    "    for i in range(len(x)):\n",
    "        A[i] = P[0] * f1(x[i]) + P[1] * f2(x[i]) + P[2] * f3(x[i])\n",
    "    return A"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using [`scipy.integrate.quad`](https://docs.scipy.org/doc/scipy/reference/generated/scipy.integrate.quad.html), we are going to check if this function is a density function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "Pf1 = quad(lambda x : P[0] * f1(x), 0, 1)\n",
    "Pf2 = quad(lambda x : P[1] * f2(x), 1, 2)\n",
    "Pf3 = quad(lambda x : P[2] * f3(x), 2, 3)\n",
    "print(f\"Integral of f1 between 0 and 1: {Pf1[0]}, with error {Pf1[1]}\")\n",
    "print(f\"Integral of f2 between 1 and 2: {Pf2[0]}, with error {Pf2[1]}\")\n",
    "print(f\"Integral of f3 between 2 and 3: {Pf3[0]}, with error {Pf3[1]}\")\n",
    "print(f\"Sum of integrals for integrating f between 0 and 3: {Pf1[0] + Pf2[0] + Pf3[0]}\")\n",
    "print(f\"Maximum of error: {max(Pf1[1], Pf2[1], Pf3[1])}\")"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "We are now going to plot this density function."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "plt.figure(figsize=(10, 7))\n",
    "plt.plot(np.linspace(0, 3, 1000), f(np.linspace(0, 3, 1000)))\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using this density function, we are now going to generate samples following this law.\n",
    "In our example, we want a sample of size $n=10000$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "V = np.random.random(10000)\n",
    "\n",
    "sample =[]\n",
    "for v in V:\n",
    "    if v <= Pf1[0]:\n",
    "        for x in np.linspace(0,1,1000):\n",
    "            if v >= quad(lambda x : P[0] * f1(x), 0, x)[0]:\n",
    "                sample.append(x)\n",
    "                break\n",
    "    elif Pf1[0] < v <= Pf1[0] + Pf2[0]:\n",
    "        for x in np.linspace(1,2,1000):\n",
    "            if v >= Pf1[0] + quad(lambda x : P[1] * f2(x), 1, x)[0]:\n",
    "                sample.append(x)\n",
    "                break\n",
    "    else:\n",
    "        for x in np.linspace(2,3,1000):\n",
    "            if v >= Pf1[0] + Pf2[0] + quad(lambda x : P[2] * f3(x), 2, x)[0]:\n",
    "                sample.append(x)\n",
    "                break\n",
    "\n",
    "plt.hist(sample, density=True, ec='white')\n",
    "plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Gaussian mixture\n",
    "\n",
    "We now consider a gaussian mixture, which is defined as follows:\n",
    "$$f(x)=p_1f_1(x)+p_2f_2(x)$$\n",
    "where $f_1$ and $f_2$ are two gaussian laws of parameters $(-3,1)$ and $(3,1)$, and $p_1$ and $p_2$ are two parameters."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "f1 = lambda x : 1 / np.sqrt(2 * np.pi) * np.exp(-(x + 3)**2/2)\n",
    "f2 = lambda x : 1 / np.sqrt(2 * np.pi) * np.exp(-(x - 3)**2/2)\n",
    "\n",
    "P = [(0.5,0.5),(1/4,3/4),(3/4,1/4)]\n",
    "\n",
    "for p in P:\n",
    "    plt.figure(figsize=(10,10))\n",
    "    plt.plot(np.linspace(-10, 10, 100), f1(np.linspace(-10, 10, 100)), color='red',label='f1')\n",
    "    plt.plot(np.linspace(-10, 10, 100), f2(np.linspace(-10, 10, 100)), color='blue', label='f2')\n",
    "    plt.plot(np.linspace(-10, 10, 100), p[0] * f1(np.linspace(-10, 10, 100)) + p[1] * f2(np.linspace(-10, 10, 100)), color='green', label='f1 + f2')\n",
    "    plt.legend()\n",
    "    plt.show()"
   ]
  },
  {
   "attachments": {},
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Using the $f$ function, we are now going to simulate a sample of size 10000 and plot their histogram.\n",
    "This will be done for each tuple of parameters $(p_1,p_2)$."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "for p in P:\n",
    "    mixed_sample = np.zeros(10000)\n",
    "\n",
    "    for i in range(10000):\n",
    "        u = np.random.random()\n",
    "        if u <= p[0]:\n",
    "            mixed_sample[i] = np.random.normal(-3, 1)\n",
    "        else:\n",
    "            mixed_sample[i] = np.random.normal(3, 1)\n",
    "\n",
    "    plt.figure(figsize=(10,10))\n",
    "    plt.hist(mixed_sample, density=True, ec='white')\n",
    "    plt.show()"
   ]
  }
 ],
 "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.9.7"
  },
  "orig_nbformat": 4,
  "vscode": {
   "interpreter": {
    "hash": "e13816212fbf99b256aa697e4560deaee36e03de1f48b7ff11501fca2851863b"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}