imp_sample.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "id": "6ecbcf28",
   "metadata": {},
   "source": [
    "# Computing Mean of a Likelihood Ratio Process"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f399a2b0",
   "metadata": {},
   "source": [
    "## Contents\n",
    "\n",
    "- [Computing Mean of a Likelihood Ratio Process](#Computing-Mean-of-a-Likelihood-Ratio-Process)  \n",
    "  - [Overview](#Overview)  \n",
    "  - [Mathematical Expectation of Likelihood Ratio](#Mathematical-Expectation-of-Likelihood-Ratio)  \n",
    "  - [Importance sampling](#Importance-sampling)  \n",
    "  - [Selecting a Sampling Distribution](#Selecting-a-Sampling-Distribution)  \n",
    "  - [Approximating a cumulative likelihood ratio](#Approximating-a-cumulative-likelihood-ratio)  \n",
    "  - [Distribution of  Sample Mean](#Distribution-of--Sample-Mean)  \n",
    "  - [More Thoughts about Choice of Sampling Distribution](#More-Thoughts-about-Choice-of-Sampling-Distribution)  "
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ba3876bc",
   "metadata": {},
   "source": [
    "## Overview\n",
    "\n",
    "In  [this lecture](https://python.quantecon.org/likelihood_ratio_process.html) we described a peculiar property of a likelihood ratio process, namely, that it’s mean equals one for all $ t \\geq 0 $ despite it’s converging to zero almost surely.\n",
    "\n",
    "While it is easy to verify that peculiar properly analytically (i.e., in population), it is challenging to use a computer simulation to verify it via an application of a law of large numbers that entails studying sample averages of repeated simulations.\n",
    "\n",
    "To confront this challenge, this lecture puts **importance sampling** to work to accelerate convergence of sample averages to population means.\n",
    "\n",
    "We use  importance sampling to estimate the mean of a cumulative likelihood ratio $ L\\left(\\omega^t\\right) = \\prod_{i=1}^t \\ell \\left(\\omega_i\\right) $.\n",
    "\n",
    "We start by importing some Python packages."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1f7fd2b1",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numba import njit, vectorize, prange\n",
    "import matplotlib.pyplot as plt\n",
    "from math import gamma"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "888827c1",
   "metadata": {},
   "source": [
    "## Mathematical Expectation of Likelihood Ratio\n",
    "\n",
    "In [this lecture](https://python.quantecon.org/likelihood_ratio_process.html), we studied a likelihood ratio $ \\ell \\left(\\omega_t\\right) $\n",
    "\n",
    "$$\n",
    "\\ell \\left( \\omega_t \\right) = \\frac{f\\left(\\omega_t\\right)}{g\\left(\\omega_t\\right)}\n",
    "$$\n",
    "\n",
    "where $ f $ and $ g $ are densities for Beta distributions with parameters $ F_a $, $ F_b $, $ G_a $, $ G_b $.\n",
    "\n",
    "Assume that an i.i.d. random variable $ \\omega_t \\in \\Omega $ is generated by $ g $.\n",
    "\n",
    "The **cumulative likelihood ratio**  $ L \\left(\\omega^t\\right) $ is\n",
    "\n",
    "$$\n",
    "L\\left(\\omega^t\\right) = \\prod_{i=1}^t \\ell \\left(\\omega_i\\right)\n",
    "$$\n",
    "\n",
    "Our goal is to approximate the mathematical expectation $ E \\left[ L\\left(\\omega^t\\right) \\right] $ well.\n",
    "\n",
    "In [this lecture](https://python.quantecon.org/likelihood_ratio_process.html), we showed that  $ E \\left[ L\\left(\\omega^t\\right) \\right] $  equals $ 1 $ for all $ t $.\n",
    "We want to check out how well this holds if we replace $ E $ by with  sample averages from simulations.\n",
    "\n",
    "This turns out to be easier said than done because for\n",
    "Beta distributions assumed above, $ L\\left(\\omega^t\\right) $ has\n",
    "a very skewed distribution with a very long tail as $ t \\rightarrow \\infty $.\n",
    "\n",
    "This property makes it difficult  efficiently and accurately to estimate the mean by standard Monte Carlo simulation methods.\n",
    "\n",
    "In this lecture we explore how a standard Monte Carlo method fails and how **importance sampling**\n",
    "provides a more computationally efficient way to approximate the mean of the cumulative likelihood ratio.\n",
    "\n",
    "We first take a look at the density functions `f` and `g` ."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "313e105d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# Parameters in the two beta distributions.\n",
    "F_a, F_b = 1, 1\n",
    "G_a, G_b = 3, 1.2\n",
    "\n",
    "@vectorize\n",
    "def p(w, a, b):\n",
    "    r = gamma(a + b) / (gamma(a) * gamma(b))\n",
    "    return r * w ** (a-1) * (1 - w) ** (b-1)\n",
    "\n",
    "# The two density functions.\n",
    "f = njit(lambda w: p(w, F_a, F_b))\n",
    "g = njit(lambda w: p(w, G_a, G_b))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "3b8d4a46",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "w_range = np.linspace(1e-5, 1-1e-5, 1000)\n",
    "\n",
    "plt.plot(w_range, g(w_range), label='g')\n",
    "plt.plot(w_range, f(w_range), label='f')\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.legend()\n",
    "plt.title('density functions $f$ and $g$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de06f460",
   "metadata": {},
   "source": [
    "The likelihood ratio is `l(w)=f(w)/g(w)`."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "904458c4",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "l = njit(lambda w: f(w) / g(w))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "56007d93",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "plt.plot(w_range, l(w_range))\n",
    "plt.title('$\\ell(\\omega)$')\n",
    "plt.xlabel('$\\omega$')\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "225b3310",
   "metadata": {},
   "source": [
    "The above plots shows that as $ \\omega \\rightarrow 0 $, $ f \\left(\\omega\\right) $ is unchanged and $ g \\left(\\omega\\right) \\rightarrow 0 $, so  the likelihood ratio approaches infinity.\n",
    "\n",
    "A Monte Carlo approximation of $ \\hat{E} \\left[L\\left(\\omega^t\\right)\\right] = \\hat{E} \\left[\\prod_{i=1}^t \\ell \\left(\\omega_i\\right)\\right] $  would repeatedly draw $ \\omega $ from $ g $, calculate the likelihood ratio $ \\ell(\\omega) = \\frac{f(\\omega)}{g(\\omega)} $ for each draw, then average these over all draws.\n",
    "\n",
    "Because $ g(\\omega) \\rightarrow 0 $ as $ \\omega \\rightarrow 0 $, such a simulation procedure  undersamples a part of the sample space $ [0,1] $ that it is important to visit often in order to do a good job of approximating the mathematical expectation of the likelihood ratio $ \\ell(\\omega) $.\n",
    "\n",
    "We illustrate this numerically below."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "13132a08",
   "metadata": {},
   "source": [
    "## Importance sampling\n",
    "\n",
    "We circumvent the issue by using a *change of distribution* called **importance sampling**.\n",
    "\n",
    "Instead of drawing from $ g $ to generate data during the simulation, we use  an alternative\n",
    "distribution $ h $ to generate draws of $ \\omega $.\n",
    "\n",
    "The idea is to design $ h $ so  that it oversamples the region of $ \\Omega $ where\n",
    "$ \\ell \\left(\\omega_t\\right) $ has large values but  low density under $ g $.\n",
    "\n",
    "After we construct a sample in this way, we must then weight each realization by the likelihood ratio of $ g $ and $ h $ when we compute the empirical mean\n",
    "of the likelihood ratio.\n",
    "\n",
    "By doing this, we properly account for the fact that we are using $ h $ and not  $ g $ to simulate data.\n",
    "\n",
    "To illustrate, suppose were interested in $ {E}\\left[\\ell\\left(\\omega\\right)\\right] $.\n",
    "\n",
    "We could simply compute:\n",
    "\n",
    "$$\n",
    "\\hat{E}^g \\left[\\ell\\left(\\omega\\right)\\right] = \\frac{1}{N} \\sum_{i=1}^{N} \\ell(w_i^g)\n",
    "$$\n",
    "\n",
    "where $ \\omega_i^g $ indicates that  $ \\omega_i $ is drawn from $ g $.\n",
    "\n",
    "But using our insight from importance sampling, we could instead calculate the object:\n",
    "\n",
    "$$\n",
    "\\hat{E}^h \\left[\\ell\\left(\\omega\\right) \\frac{g(w)}{h(w)} \\right] = \\frac{1}{N} \\sum_{i=1}^{N} \\ell(w_i^h) \\frac{g(w_i^h)}{h(w_i^h)}\n",
    "$$\n",
    "\n",
    "where $ w_i $ is now drawn from importance distribution $ h $.\n",
    "\n",
    "Notice that the above two  are exactly the same population objects:\n",
    "\n",
    "$$\n",
    "E^g\\left[\\ell\\left(\\omega\\right)\\right] = \\int_\\Omega \\ell(\\omega) g(\\omega) d\\omega = \\int_\\Omega \\ell(\\omega) \\frac{g(\\omega)}{h(\\omega)} h(\\omega) d\\omega = E^h\\left[\\ell\\left(\\omega\\right) \\frac{g(\\omega)}{h(\\omega)}\\right]\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "57d421a0",
   "metadata": {},
   "source": [
    "## Selecting a Sampling Distribution\n",
    "\n",
    "Since we must use an $ h $ that has larger mass in parts of the distribution to which  $ g $ puts low mass, we use $ h=Beta(0.5, 0.5) $ as our importance distribution.\n",
    "\n",
    "The  plots compare $ g $ and $ h $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "39b14546",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "g_a, g_b = G_a, G_b\n",
    "h_a, h_b = 0.5, 0.5"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ec1652d2",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "w_range = np.linspace(1e-5, 1-1e-5, 1000)\n",
    "\n",
    "plt.plot(w_range, g(w_range), label=f'g=Beta({g_a}, {g_b})')\n",
    "plt.plot(w_range, p(w_range, 0.5, 0.5), label=f'h=Beta({h_a}, {h_b})')\n",
    "plt.title('real data generating process $g$ and importance distribution $h$')\n",
    "plt.legend()\n",
    "plt.ylim([0., 3.])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a5c5a1d8",
   "metadata": {},
   "source": [
    "## Approximating a cumulative likelihood ratio\n",
    "\n",
    "We now study how to use importance sampling to approximate\n",
    "$ {E} \\left[L(\\omega^t)\\right] = \\left[\\prod_{i=1}^T \\ell \\left(\\omega_i\\right)\\right] $.\n",
    "\n",
    "As above, our plan is to draw  sequences $ \\omega^t $ from $ q $ and then re-weight the likelihood ratio appropriately:\n",
    "\n",
    "$$\n",
    "\\hat{E}^p \\left[L\\left(\\omega^t\\right)\\right] = \\hat{E}^p \\left[\\prod_{t=1}^T \\ell \\left(\\omega_t\\right)\\right] = \\hat{E}^q \\left[\\prod_{t=1}^T \\ell \\left(\\omega_t\\right) \\frac{p\\left(\\omega_{t}\\right)}{q\\left(\\omega_{t}\\right)}\\right] =\n",
    "\\frac{1}{N} \\sum_{i=1}^{N}\\left( \\prod_{t=1}^{T} \\ell(\\omega_{i,t}^h)\\frac{p\\left(\\omega_{i,t}^h\\right)}{q\\left(\\omega_{i,t}^h\\right)}\\right)\n",
    "$$\n",
    "\n",
    "where the last equality uses $ \\omega_{i,t}^h $ drawn from the importance distribution $ q $.\n",
    "\n",
    "Here $ \\frac{p\\left(\\omega_{i,t}^q\\right)}{q\\left(\\omega_{i,t}^q\\right)} $ is the weight we assign to each data point $ \\omega_{i,t}^q $.\n",
    "\n",
    "Below we prepare a Python function for computing the importance sampling estimates given any beta distributions $ p $, $ q $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "05400d77",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit(parallel=True)\n",
    "def estimate(p_a, p_b, q_a, q_b, T=1, N=10000):\n",
    "\n",
    "    μ_L = 0\n",
    "    for i in prange(N):\n",
    "\n",
    "        L = 1\n",
    "        weight = 1\n",
    "        for t in range(T):\n",
    "            w = np.random.beta(q_a, q_b)\n",
    "            l = f(w) / g(w)\n",
    "\n",
    "            L *= l\n",
    "            weight *= p(w, p_a, p_b) / p(w, q_a, q_b)\n",
    "\n",
    "        μ_L += L * weight\n",
    "\n",
    "    μ_L /= N\n",
    "\n",
    "    return μ_L"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "09d49b64",
   "metadata": {},
   "source": [
    "Consider the case when $ T=1 $, which amounts  to approximating $ E_0\\left[\\ell\\left(\\omega\\right)\\right] $\n",
    "\n",
    "For the standard Monte Carlo estimate, we can set $ p=g $ and $ q=g $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "646e58bc",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "estimate(g_a, g_b, g_a, g_b, T=1, N=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "9a1d22b9",
   "metadata": {},
   "source": [
    "For our importance sampling estimate, we set $ q = h $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "26cb0e19",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "estimate(g_a, g_b, h_a, h_b, T=1, N=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4390ef57",
   "metadata": {},
   "source": [
    "Evidently, even at T=1, our importance sampling  estimate is closer to $ 1 $ than is the Monte Carlo estimate.\n",
    "\n",
    "Bigger differences arise when computing expectations over longer sequences, $ E_0\\left[L\\left(\\omega^t\\right)\\right] $.\n",
    "\n",
    "Setting $ T=10 $, we find that the  Monte Carlo method severely underestimates the mean while importance sampling\n",
    "still produces an estimate close to its theoretical value of unity."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "6774ab5d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "estimate(g_a, g_b, g_a, g_b, T=10, N=10000)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d8caf41d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "estimate(g_a, g_b, h_a, h_b, T=10, N=10000)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "605efbe4",
   "metadata": {},
   "source": [
    "## Distribution of  Sample Mean\n",
    "\n",
    "We next study the bias and efficiency of the Monte Carlo and importance sampling approaches.\n",
    "\n",
    "The code  below produces distributions of estimates using both Monte Carlo and importance sampling methods."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "d76f2c1d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "@njit(parallel=True)\n",
    "def simulate(p_a, p_b, q_a, q_b, N_simu, T=1):\n",
    "\n",
    "    μ_L_p = np.empty(N_simu)\n",
    "    μ_L_q = np.empty(N_simu)\n",
    "\n",
    "    for i in prange(N_simu):\n",
    "        μ_L_p[i] = estimate(p_a, p_b, p_a, p_b, T=T)\n",
    "        μ_L_q[i] = estimate(p_a, p_b, q_a, q_b, T=T)\n",
    "\n",
    "    return μ_L_p, μ_L_q"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f1871492",
   "metadata": {},
   "source": [
    "Again, we first consider estimating $ {E} \\left[\\ell\\left(\\omega\\right)\\right] $ by setting T=1.\n",
    "\n",
    "We simulate $ 1000 $ times for each method."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "2b138498",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "N_simu = 1000\n",
    "μ_L_p, μ_L_q = simulate(g_a, g_b, h_a, h_b, N_simu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "1556b285",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# standard Monte Carlo (mean and std)\n",
    "np.nanmean(μ_L_p), np.nanvar(μ_L_p)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "de778b6d",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# importance sampling (mean and std)\n",
    "np.nanmean(μ_L_q), np.nanvar(μ_L_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "aad650f7",
   "metadata": {},
   "source": [
    "Although both methods tend to provide a mean estimate of $ {E} \\left[\\ell\\left(\\omega\\right)\\right] $ close to $ 1 $, the importance sampling estimates have  smaller variance.\n",
    "\n",
    "Next, we present distributions of estimates for $ \\hat{E} \\left[L\\left(\\omega^t\\right)\\right] $, in cases for $ T=1, 5, 10, 20 $."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5c00e266",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "fig, axs = plt.subplots(2, 2, figsize=(14, 10))\n",
    "\n",
    "μ_range = np.linspace(0, 2, 100)\n",
    "\n",
    "for i, t in enumerate([1, 5, 10, 20]):\n",
    "    row = i // 2\n",
    "    col = i % 2\n",
    "\n",
    "    μ_L_p, μ_L_q = simulate(g_a, g_b, h_a, h_b, N_simu, T=t)\n",
    "    μ_hat_p, μ_hat_q = np.nanmean(μ_L_p), np.nanmean(μ_L_q)\n",
    "    σ_hat_p, σ_hat_q = np.nanvar(μ_L_p), np.nanvar(μ_L_q)\n",
    "\n",
    "    axs[row, col].set_xlabel('$μ_L$')\n",
    "    axs[row, col].set_ylabel('frequency')\n",
    "    axs[row, col].set_title(f'$T$={t}')\n",
    "    n_p, bins_p, _ = axs[row, col].hist(μ_L_p, bins=μ_range, color='r', alpha=0.5, label='$g$ generating')\n",
    "    n_q, bins_q, _ = axs[row, col].hist(μ_L_q, bins=μ_range, color='b', alpha=0.5, label='$h$ generating')\n",
    "    axs[row, col].legend(loc=4)\n",
    "\n",
    "    for n, bins, μ_hat, σ_hat in [[n_p, bins_p, μ_hat_p, σ_hat_p],\n",
    "                                  [n_q, bins_q, μ_hat_q, σ_hat_q]]:\n",
    "        idx = np.argmax(n)\n",
    "        axs[row, col].text(bins[idx], n[idx], '$\\hat{μ}$='+f'{μ_hat:.4g}'+', $\\hat{σ}=$'+f'{σ_hat:.4g}')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "27512a39",
   "metadata": {},
   "source": [
    "The simulation exercises above show that the importance sampling estimates are unbiased under all $ T $\n",
    "while the standard Monte Carlo estimates are biased downwards.\n",
    "\n",
    "Evidently, the bias increases with increases in $ T $."
   ]
  },
  {
   "cell_type": "markdown",
   "id": "542c8561",
   "metadata": {},
   "source": [
    "## More Thoughts about Choice of Sampling Distribution\n",
    "\n",
    "Above, we arbitraily chose $ h = Beta(0.5,0.5) $ as the importance distribution.\n",
    "\n",
    "Is there an optimal importance distribution?\n",
    "\n",
    "In our particular case, since we  know in advance that $ E_0 \\left[ L\\left(\\omega^t\\right) \\right] = 1 $.\n",
    "\n",
    "We can use that knowledge to our advantage.\n",
    "\n",
    "Thus, suppose that we simply use  $ h = f $.\n",
    "\n",
    "When estimating the mean of the likelihood ratio (T=1), we get:\n",
    "\n",
    "$$\n",
    "\\hat{E}^f \\left[\\ell(\\omega) \\frac{g(\\omega)}{f(\\omega)} \\right] = \\hat{E}^f \\left[\\frac{f(\\omega)}{g(\\omega)} \\frac{g(\\omega)}{f(\\omega)} \\right] = \\frac{1}{N} \\sum_{i=1}^{N} \\ell(w_i^f) \\frac{g(w_i^f)}{f(w_i^f)} = 1\n",
    "$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "119a3d09",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "μ_L_p, μ_L_q = simulate(g_a, g_b, F_a, F_b, N_simu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "767662bf",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "# importance sampling (mean and std)\n",
    "np.nanmean(μ_L_q), np.nanvar(μ_L_q)"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "59d88a5c",
   "metadata": {},
   "source": [
    "We could also use other distributions as our importance distribution.\n",
    "\n",
    "Below we choose just a few and compare their sampling properties."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "79fa3884",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "a_list = [0.5, 1., 2.]\n",
    "b_list = [0.5, 1.2, 5.]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "b10baa79",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "w_range = np.linspace(1e-5, 1-1e-5, 1000)\n",
    "\n",
    "plt.plot(w_range, g(w_range), label=f'p=Beta({g_a}, {g_b})')\n",
    "plt.plot(w_range, p(w_range, a_list[0], b_list[0]), label=f'g=Beta({a_list[0]}, {b_list[0]})')\n",
    "plt.plot(w_range, p(w_range, a_list[1], b_list[1]), label=f'g=Beta({a_list[1]}, {b_list[1]})')\n",
    "plt.plot(w_range, p(w_range, a_list[2], b_list[2]), label=f'g=Beta({a_list[2]}, {b_list[2]})')\n",
    "plt.title('real data generating process $g$ and importance distribution $h$')\n",
    "plt.legend()\n",
    "plt.ylim([0., 3.])\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "4c766137",
   "metadata": {},
   "source": [
    "We consider two additonal distributions.\n",
    "\n",
    "As a reminder $ h_1 $ is the original $ Beta(0.5,0.5) $ distribution that we used above.\n",
    "\n",
    "$ h_2 $ is the $ Beta(1,1.2) $ distribution.\n",
    "\n",
    "Note how $ h_2 $ has a similar shape to $ g $ at higher values of distribution but more mass at lower values.\n",
    "\n",
    "Our hunch is that $ h_2 $ should be a good importance sampling distribution.\n",
    "\n",
    "$ h_3 $ is the $ Beta(2,5) $ distribution.\n",
    "\n",
    "Note how $ h_3 $ has zero mass at values very close to 0 and at values close to 1.\n",
    "\n",
    "Our hunch is that $ h_3 $ will  be a poor importance sampling distribution.\n",
    "\n",
    "We first simulate a plot the distribution of estimates for $ \\hat{E} \\left[L\\left(\\omega^t\\right)\\right] $ using $ h_2 $ as the importance sampling distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "ca60f018",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "h_a = a_list[1]\n",
    "h_b = b_list[1]\n",
    "\n",
    "fig, axs = plt.subplots(1,2, figsize=(14, 10))\n",
    "\n",
    "μ_range = np.linspace(0, 2, 100)\n",
    "\n",
    "for i, t in enumerate([1, 20]):\n",
    "\n",
    "\n",
    "    μ_L_p, μ_L_q = simulate(g_a, g_b, h_a, h_b, N_simu, T=t)\n",
    "    μ_hat_p, μ_hat_q = np.nanmean(μ_L_p), np.nanmean(μ_L_q)\n",
    "    σ_hat_p, σ_hat_q = np.nanvar(μ_L_p), np.nanvar(μ_L_q)\n",
    "\n",
    "    axs[i].set_xlabel('$μ_L$')\n",
    "    axs[i].set_ylabel('frequency')\n",
    "    axs[i].set_title(f'$T$={t}')\n",
    "    n_p, bins_p, _ = axs[i].hist(μ_L_p, bins=μ_range, color='r', alpha=0.5, label='$g$ generating')\n",
    "    n_q, bins_q, _ = axs[i].hist(μ_L_q, bins=μ_range, color='b', alpha=0.5, label='$h_2$ generating')\n",
    "    axs[i].legend(loc=4)\n",
    "\n",
    "    for n, bins, μ_hat, σ_hat in [[n_p, bins_p, μ_hat_p, σ_hat_p],\n",
    "                                  [n_q, bins_q, μ_hat_q, σ_hat_q]]:\n",
    "        idx = np.argmax(n)\n",
    "        axs[i].text(bins[idx], n[idx], '$\\hat{μ}$='+f'{μ_hat:.4g}'+', $\\hat{σ}=$'+f'{σ_hat:.4g}')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "2af359fb",
   "metadata": {},
   "source": [
    "Our simulations suggest that indeed $ h_2 $ is a quite  good importance sampling distribution for our problem.\n",
    "\n",
    "Even at $ T=20 $, the mean  is very close to $ 1 $ and the  variance is small."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "5e7dc944",
   "metadata": {
    "hide-output": false
   },
   "outputs": [],
   "source": [
    "h_a = a_list[2]\n",
    "h_b = b_list[2]\n",
    "\n",
    "fig, axs = plt.subplots(1,2, figsize=(14, 10))\n",
    "\n",
    "μ_range = np.linspace(0, 2, 100)\n",
    "\n",
    "for i, t in enumerate([1, 20]):\n",
    "\n",
    "\n",
    "    μ_L_p, μ_L_q = simulate(g_a, g_b, h_a, h_b, N_simu, T=t)\n",
    "    μ_hat_p, μ_hat_q = np.nanmean(μ_L_p), np.nanmean(μ_L_q)\n",
    "    σ_hat_p, σ_hat_q = np.nanvar(μ_L_p), np.nanvar(μ_L_q)\n",
    "\n",
    "    axs[i].set_xlabel('$μ_L$')\n",
    "    axs[i].set_ylabel('frequency')\n",
    "    axs[i].set_title(f'$T$={t}')\n",
    "    n_p, bins_p, _ = axs[i].hist(μ_L_p, bins=μ_range, color='r', alpha=0.5, label='$g$ generating')\n",
    "    n_q, bins_q, _ = axs[i].hist(μ_L_q, bins=μ_range, color='b', alpha=0.5, label='$h_3$ generating')\n",
    "    axs[i].legend(loc=4)\n",
    "\n",
    "    for n, bins, μ_hat, σ_hat in [[n_p, bins_p, μ_hat_p, σ_hat_p],\n",
    "                                  [n_q, bins_q, μ_hat_q, σ_hat_q]]:\n",
    "        idx = np.argmax(n)\n",
    "        axs[i].text(bins[idx], n[idx], '$\\hat{μ}$='+f'{μ_hat:.4g}'+', $\\hat{σ}=$'+f'{σ_hat:.4g}')\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7602b9e8",
   "metadata": {},
   "source": [
    "However, $ h_3 $ is evidently a poor importance sampling distribution forpir problem,\n",
    "with a mean estimate far away from $ 1 $ for $ T = 20 $.\n",
    "\n",
    "Notice that evan at $ T = 1 $, the mean estimate with importance sampling is more biased than just sampling with  $ g $ itself.\n",
    "\n",
    "Thus, our simulations suggest that we would be better off simply using Monte Carlo\n",
    "approximations under $ g $ than using $ h_3 $ as an importance sampling distribution for our problem."
   ]
  }
 ],
 "metadata": {
  "date": 1714442505.78251,
  "filename": "imp_sample.md",
  "kernelspec": {
   "display_name": "Python",
   "language": "python3",
   "name": "python3"
  },
  "title": "Computing Mean of a Likelihood Ratio Process"
 },
 "nbformat": 4,
 "nbformat_minor": 5
}