模拟精度改进技术.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "<h1> <center> 模拟精度改进技术 </h1>\n",
    "\n",
    "# 理论\n",
    "\n",
    "## 引言\n",
    "\n",
    "- 在模拟研究中, 人们感兴趣的是想确定某一个参数$\\theta$(来自于某一个随机模型)\n",
    "- 为估计$\\theta$, 通过对模型进行模拟得到输出的数据$X$, 满足$\\theta=E[X]$\n",
    "- 重复模拟, 直到n次后停止, 记第$i$次的输出的变量为$X_i$, 此时$\\theta$的估计为$\\bar X = \\sum_{i=1}^n \\frac{X_i}{n}$\n",
    "- 该估计无偏, 则均方误差为其方差\n",
    "\n",
    "$$\n",
    "\\text{MSE} = E[(\\bar X -\\theta)^2] = \\text{Var}(\\bar{X}) = \\text{Var}(X)/n\n",
    "$$\n",
    "\n",
    "- 若能得到比$\\bar{X}$的方差还要小的 $\\theta$的无偏估计, 就可以对$\\theta$的估计进行改进\n",
    "\n",
    "## 对偶变量法\n",
    "\n",
    "假设我们对估计$\\theta = E[X]$感兴趣, 有假设产生了同分布且具有均值$\\theta$的随机变量$X_1, X_2$, 则:\n",
    "\n",
    "$$\n",
    "\\operatorname{Var}\\left(\\frac{X_{1}+X_{2}}{2}\\right)=\\frac{1}{4}\\left[\\operatorname{Var}\\left(X_{1}\\right)+\\operatorname{Var}\\left(X_{2}\\right)+2 \\operatorname{Cov}\\left(X_{1}, X_{2}\\right)\\right]\n",
    "$$\n",
    "\n",
    "此时易知若$X_1,X_2$负相关, 则有$\\operatorname{Var}\\left(\\frac{X_{1}+X_{2}}{2}\\right) \\leqslant \\frac{1}{4}\\left[\\operatorname{Var}\\left(X_{1}\\right)+\\operatorname{Var}\\left(X_{2}\\right)\\right]$, 则有在某种意义上来说方差减小了, 也就是, 估计的精度得到提高, 而仍然有$E[\\frac{X_1+X_2}{2}] = \\theta$.\n",
    "\n",
    "以下<font color=red>**定理**</font>可以使得我们的结论成立, 也就是说我们只要保证两变量为负相关就可以找到相关的结果. 此处为找对偶变量的一种方法.\n",
    "\n",
    "> 如果 $h\\left(x_{1}, x_{2}, \\cdots, x_{n}\\right)$ 是其每个自变量的单调函数, 则对独立随机数的集合 $U_{1}, U_{2}, \\cdots, U_{n}$, 有 \n",
    "> $$\\operatorname{Cov}\\left[h\\left(U_{1}, U_{2}, \\cdots, U_{n}\\right), h\\left(1-U_{1}, 1-U_{2}, \\cdots, 1-U_{n}\\right)\\right] \\leqslant 0$$\n",
    "\n",
    "---\n",
    "\n",
    "### 例 6.1\n",
    "\n",
    "假设需要估计$\\theta = \\int_0^1 e^x dx = e-1$. 显然 $\\theta=E[e^U]$, 此时有:\n",
    "\n",
    "$$\n",
    "\\operatorname{Var}(e^U) = E[e^{2U}]-E[e^U] ^2 = \\frac{e^2-1}{2} - (e-1)^2\n",
    "$$\n",
    "\n",
    "- 使用独立抽样法:\n",
    "  - 产生两个随机数$U_1,U_2$, 令$\\frac{e^{U_1}+e^{U_2}}{2}$为$\\theta$的估计. 此时该估计为无偏估计\n",
    "  - 方差为: $\\operatorname{Var}(\\frac{e^{U_1}+e^{U_2}}{2})=\\frac{\\operatorname{Var}(e^U)}{2} \\approx 0.1210$\n",
    "- 使用对偶变量法:\n",
    "  - 此时首先要论证对偶变量法从理论上优于独立变量法的. \n",
    "    - 此时的`h`是一个指数函数, 得出单调. 则从理论上说明应该是优于独立抽样法的\n",
    "  - 计算: $\\operatorname{Var}(\\frac{e^U+e^{1-U}}{2})=\\frac{\\operatorname{Var}(e^U)}{2}+\\frac{\\operatorname{Cov}(e^U, e^{1-U})}{2}\\approx 0.0039$\n",
    "\n",
    "综上, 结果变好了\n",
    "\n",
    "---\n",
    "\n",
    "### 例 6.2\n",
    "\n",
    "令$X_j, j=1, 2, \\cdots, 5$是均值为1的独立的指数随机变量, 模拟估计:\n",
    "\n",
    "$$\n",
    "\\theta = P\\{\n",
    "    \\sum_{j=1}^5 j X_j \\ge 21.6\n",
    "    \\}\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "若用对偶变量法对 $\\theta$ 进行估计, 我们按如下方式产生 $2 n$ 组随机数:\n",
    "- (奇数列) 独立产生 $n$ 组随机数: $U_{1}^{2 k-1}, U_{2}^{2 k-1}, \\cdots, U_{k}^{2 k-1}, k=1,2, \\cdots, n$\n",
    "- (偶数列) 令 $U_{j}^{2 k}=1-U_{j}^{2 k-1}, j=1,2, \\cdots, 5, k=1,2, \\cdots, n$\n",
    "令 $X_{j}^{i}=-\\ln \\left(U_{j}^{i}\\right), j=1,2, \\cdots, 5 ; i=1,2, \\cdots, 2 n$\n",
    "\n",
    "$$\n",
    "\\tilde{I}_{i}=\\left\\{\\begin{array}{cc}\n",
    "1, & \\text { 若 } \\sum_{j=1}^{5} j X_{j}^{i} \\geq 21.6 \\\\\n",
    "0, & \\text { 否则. }\n",
    "\\end{array}\\right.\n",
    "$$\n",
    "\n",
    "$\\left\\{\\tilde{I}_{i}\\right\\}$ 奇偶列各自成员间相互独立, 但相邻奇偶数不一定独立\n",
    "- $\\left\\{\\tilde{I}_{1}, \\tilde{I}_{3}, \\ldots\\right\\}$ 独立; $\\left\\{\\tilde{I}_{2}, \\tilde{I}_{4}, \\ldots\\right\\}$ 独立; $\\left\\{\\left(\\tilde{I}_{1}, \\tilde{I}_{2}\\right),\\left(\\tilde{I}_{3}, \\tilde{I}_{4}\\right), \\ldots\\right\\}$ 独立; $\\tilde{I}_{1}$ 和 $\\tilde{I}_{2}$ 不独立, ...\n",
    "用 $\\bar{I}_{2}=\\sum_{i=1}^{2 n} \\frac{\\tilde{I}_{i}}{2 n}$ 对 $\\theta$ 进行估计. (频率估计概率原理)\n",
    "\n",
    "容易证明 $E\\left(\\tilde{I}_{i}\\right)=\\theta, \\operatorname{Var}\\left(\\tilde{I}_{i}\\right)=\\theta(1-\\theta), E\\left(\\bar{I}_{2}\\right)=\\theta$,\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "\\operatorname{Var}\\left(\\bar{I}_{2}\\right) &=\\frac{\\theta(1-\\theta)}{2 n}+\\frac{\\sum_{i=1}^{n} 2 \\operatorname{Cov}\\left(\\tilde{I}_{2 i}, \\tilde{I}_{2 i-1}\\right)}{4 n^{2}}=\\frac{\\theta(1-\\theta)}{2 n}+\\frac{\\operatorname{Cov}\\left(\\tilde{I}_{1} \\tilde{I}_{2}\\right)}{2 n} \\\\\n",
    "&=\\frac{\\theta(1-\\theta)}{2 n}+\\frac{E\\left(\\tilde{I}_{1} \\tilde{I}_{2}\\right)-\\theta^{2}}{2 n}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "容易看到 $I_{i}$ 是其坐标的单调函数, 得 $\\operatorname{Cov}\\left(\\tilde{I}_{1} \\tilde{I}_{2}\\right) \\leq 0$, 故 $\\bar{I}_{2}$ 的方差不大于 $\\bar{I}_{1}$ 的方差."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "from numpy import random\n",
    "from math import *\n",
    "from scipy import integrate"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# 独立抽样法\n",
    "def eg6_2_1(size = 10000):\n",
    "    # 使其生成的随机数与对偶问题保持一致\n",
    "    size *= 2\n",
    "\n",
    "    X = np.zeros((size, 5)) # store the random number\n",
    "    \n",
    "    for i in range(size):\n",
    "        for j in range(5):\n",
    "            X[i, j] = -log(random.uniform())*(j+1) # 注意此处为k+1\n",
    "        \n",
    "    X = np.sum(X, axis=1)\n",
    "    k = np.count_nonzero(X >= 21.6)\n",
    "    X = X>=21.6\n",
    "\n",
    "    return k/size, k/size*(1-k/size)/size, np.var(X)/size\n",
    "\n",
    "# 对偶变量法\n",
    "def eg6_2_2(size = 10000):\n",
    "    X = np.zeros((size*2, 5))\n",
    "\n",
    "    for i in range(size):\n",
    "        for j in range(5):\n",
    "            X[i, j] = -log(random.uniform())\n",
    "\n",
    "    X[size:2*size,:] = -np.log(1-np.exp(-X[:size, :])) # 参见习题中的变换\n",
    "\n",
    "    X = X * np.arange(1, 6)\n",
    "\n",
    "    X = np.sum(X, axis=1)\n",
    "\n",
    "    X = X>=21.6\n",
    "    # NOTE: 此处不要出现格式问题\n",
    "    x1 = np.array(X[:size], dtype=np.float)\n",
    "    x2 = np.array(X[size:], dtype=np.float)\n",
    "    \n",
    "    r = np.count_nonzero(x1*x2)\n",
    "\n",
    "    X = (x1+x2)/2\n",
    "    th = np.mean(X)\n",
    "\n",
    "    r = float(r)\n",
    "    size = float(size)\n",
    "    \n",
    "    # NOTE: 为了验证方差的计算值的正确性, 同时计算理论值与实际计算值\n",
    "    # NOTE: 此处采用极大似然估计方差, 所以和R中的计算略有差异, 但不影响最终的结论\n",
    "    return th, np.var(X)/(size), th*(1-th)/(2*size) + (r/size-(th)**2)/(2*size)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "(0.17435, 7.197603875e-06, 7.197603875e-06)\n",
      "(0.1717, 5.776911e-06, 5.776911000000001e-06)\n"
     ]
    }
   ],
   "source": [
    "random.seed(0)\n",
    "print(eg6_2_1())\n",
    "\n",
    "random.seed(0)\n",
    "print(eg6_2_2())"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "从上述内容我们可以发现减小了方差\n",
    "\n",
    "----"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "### 对偶变量法(正态分布情形)\n",
    "\n",
    "- 提出正态分布的motivation是什么: 此时的`h`并没有一个\"显式\"的函数, 所以对于此类<font color=red>**不方便使用逆变换法的**</font>, 例如正态分布, 就可以使用以下的方式进行验证\n",
    "\n",
    "如果运用正态随机变量$Y \\sim N\\left(\\mu, \\sigma^{2}\\right)$ 模拟估计 $E\\left[h\\left(Y_{1}, Y_{2}, \\cdots, Y_{n}\\right)\\right](h$**是其坐标的单调函数时**), 我们也可以用对偶变量法减少方差.\n",
    "\n",
    "首先产生均值为 $\\mu$ 方差为 $\\sigma^{2}$ 的正态随机变量 $Y$, 然后取对偶变量 $2 \\mu-Y$, 易知\n",
    "$2 \\mu-Y$ 也具有均值 $\\mu$ 和方差 $\\sigma^{2}$. 显然它和 $Y$ 是负相关的.\n",
    "\n",
    "- $E(2 \\mu-Y)=2 \\mu-E(Y)=\\mu$\n",
    "- $\\operatorname{Var}(2 \\mu-Y)=\\operatorname{Var}(Y)=\\sigma^{2}$\n",
    "- $\\operatorname{Cov}(2 \\mu-Y, Y)=\\operatorname{Cov}(2 \\mu, Y)-\\operatorname{Cov}(Y, Y)=0-\\sigma^{2}<0$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 条件期望法\n",
    "\n",
    "### 条件期望回顾\n",
    "\n",
    "该部分直接以连续型随机变量为例进行推导(主要参考自统计推断)\n",
    "\n",
    "- 条件分布定义: 设$X, Y$是连续的二维随机变量, 其联合概率密度函数为$f(x,y)$, 边缘密度函数分别为$f_X(x), f_Y(y)$. 对任意满足$f_X(x)> 0$的$x, Y$在条件$X=x$下的条件概率密度函数(conditional pdf of $Y$ given that $X=x$)**是$y$的一个函数**, 记作$f(y|x)$, 定义为:\n",
    "  $$\n",
    "    f(y \\mid x) = \\frac{f(x, y)}{f_X(x)}\n",
    "  $$\n",
    "- 设$g(Y)$是$Y$的函数, 则$g(Y)$在$X=x$条件下的条件期望(conditional expected value of $g(Y)$ given that $X=x$)记作, 定义为:\n",
    "  $$\n",
    "    E(g(Y) \\mid x) = \\int_{-\\infty}^{+\\infty} g(y)f(y\\mid x)dy\n",
    "  $$\n",
    "\n",
    "例题:\n",
    "\n",
    "- $(X,Y)$为连续随机向量\n",
    "- $f(x,y) = e^{-y}, 0<x<y<+\\infty.$\n",
    "- 求$f(y|x)$.\n",
    "\n",
    "$$\n",
    "\\begin{aligned}\n",
    "    f_X(x) = \\int_{-\\infty}^{+\\infty} f(x, y)dy = \\int_x^{+\\infty} e^{-y}dy = e^{-x}\n",
    "\\end{aligned}\n",
    "$$\n",
    "\n",
    "$$\n",
    "\\Rightarrow \n",
    "f(y|x)\n",
    "\\begin{cases}\n",
    "    e^{-(y-x)}, &y > x > 0\\\\\n",
    "    0, &y\\le x\n",
    "\\end{cases}\n",
    "$$\n",
    "\n",
    "$$\n",
    "E(Y|X=x) = \\int_x^{+\\infty}y e^{-(y-x)}dy = 1+x\n",
    "$$\n",
    "\n",
    "根据方差公式, 可以得出:$\\text{Var}(Y|X=x) = 1$, 将Y的条件方差与无条件方差作比较, $Y\\sim \\text{Gamma}(2, 1)$, $\\text{Var}(Y) = 2$, 方差缩小.\n",
    "\n",
    "### 方差缩减原理\n",
    "\n",
    "由方差非负以及下述式子可知:\n",
    "\n",
    "$$\n",
    "\\text{Var}(X) = E[\\text{Var}(X|Y)]+\\text{Var}(E[X|Y]) \\Rightarrow \\text{Var}(X) \\ge \\text{Var}(E[X|Y])\n",
    "$$"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "\n",
    "## 重要抽样法\n",
    "\n",
    "> 写在前面: 重要抽样法主要参考[李东风老师的统计计算](https://www.math.pku.edu.cn/teachers/lidf/docs/statcomp/html/_statcompbook/sim-impsamp.html#sim-impsamp-std), 李东风老师的书写得非常详细, 此处的笔记主要介绍个人对该方法思想的理解. 重要抽样法同时在强化学习中也得到很多应用.\n",
    "\n",
    "重要抽样法主要分为两种情形:\n",
    "\n",
    "$$\n",
    "\\begin{align}\n",
    "    &I = \\int_C h(\\boldsymbol x) \\,d \\boldsymbol x \\\\\n",
    "    &E h(\\boldsymbol Y) = \\int h(\\boldsymbol y) f(\\boldsymbol y) \\,d\\boldsymbol y\n",
    "\\end{align}\n",
    "$$\n",
    "\n",
    "分别是直接求积分和求随机变量($\\boldsymbol{Y}$)的期望. 其中$\\boldsymbol Y \\sim f(\\boldsymbol y)$.\n",
    "\n",
    "- 直接求积分(红色部分非常重要, 也是对该方法进行理解的前提)\n",
    "    - 首先设<font color=red>$g(\\boldsymbol x), \\boldsymbol x \\in C$是一个概率密度的表达式, 并且其形状和$h(\\boldsymbol{x})$相近(更规范的定义可以参见其书籍主页)</font>. $X_i \\mathop{\\sim}\\limits^\\text{i.i.d} g(\\boldsymbol{x}), i=1, 2, \\cdots, N$, 假设:\n",
    "    $$\n",
    "    \\begin{aligned}\n",
    "    \\eta_i = \\frac{h(\\boldsymbol X_i)}{g(\\boldsymbol X_i)},\n",
    "    \\ i=1,2,\\dots,N\n",
    "    \\end{aligned}\n",
    "    $$\n",
    "    - 则有:\n",
    "    $$\n",
    "    \\begin{aligned}\n",
    "    E \\eta_1 = \\int_C \\frac{h(\\boldsymbol x)}{g(\\boldsymbol x)} g(\\boldsymbol x) \\,d \\boldsymbol x \n",
    "    = \\int_C h(\\boldsymbol x) \\,d \\boldsymbol x = I\n",
    "    \\end{aligned}\n",
    "    $$\n",
    "    - 此时可以使用平均值进行估计:\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    \\hat I_3 = \\bar\\eta = \\frac{1}{N} \\sum_{i=1}^N \\frac{h(\\boldsymbol X_i)}{g(\\boldsymbol X_i)}.\n",
    "    \\end{align}\n",
    "    $$\n",
    "    - 其方差为:\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    \\text{Var}(\\hat I_3) = \\text{Var}\\left( \\frac{h(\\boldsymbol X)}{g(\\boldsymbol X)} \\right) \\frac{1}{N}\n",
    "    = O(\\frac{1}{N}),\n",
    "    \\end{align}\n",
    "    $$\n",
    "    - 此时论证其方差较小的思路为:\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    \\text{Var}\\left( \\frac{h(\\boldsymbol X)}{g(\\boldsymbol X)} \\right) =& E \\left( \\frac{h^2(\\boldsymbol X)}{g^2(\\boldsymbol X)} \\right) - \\left( E \\frac{h(\\boldsymbol X)}{g(\\boldsymbol X)} \\right)^2 \n",
    "    = E \\left( \\frac{h^2(\\boldsymbol X)}{g^2(\\boldsymbol X)} \\right) - \\left( \\int_C h(\\boldsymbol x) \\,d\\boldsymbol x \\right)^2 \\\\\n",
    "    \\geq& \\left( E \\frac{|h(\\boldsymbol X)|}{g(\\boldsymbol X)} \\right)^2 - I^2\n",
    "        \\quad \\text{(Jensen's inequality)} \\\\\n",
    "    =& \\left( \\int_C |h(\\boldsymbol x)| \\,d\\boldsymbol x \\right)^2 - I^2 , \n",
    "    \\end{align}\n",
    "    $$\n",
    "    - 此时我们可以通过Jensen's inequality(取等号的情形为$\\frac{|h(\\boldsymbol X)|}{g(\\boldsymbol X)}$为一常数)可以发现我们从理论上可以将其方差缩的非常小.该部分在武大版统计计算中直接将`h`定义为＞0的函数, 从理解上来说更加方便了.\n",
    "\n",
    "- 求随机变量的期望\n",
    "    - 该部分的论证和上述比较类似, 同样定义一个试投密度$g(\\boldsymbol{x})$, 则可以得到一个无偏估计:\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    \\hat I_{3.1} =& \\frac{1}{N} \\sum_{i=1}^N h(\\boldsymbol X_i) \\frac{f(\\boldsymbol X_i)}{g(\\boldsymbol X_i)}\n",
    "    \\end{align}\n",
    "    $$\n",
    "    - 此时可以换一种写法对其进行理解, 假设$W_i = \\frac{f(\\boldsymbol X_i)}{g(\\boldsymbol X_i)}$($W_i$被称为重要性权重):\n",
    "    $$\n",
    "    \\begin{align}\n",
    "    \\hat I_{3.1} =& \\frac{1}{N} \\sum_{i=1}^N W_i h(\\boldsymbol X_i) .\n",
    "    \\end{align}\n",
    "    $$\n",
    "    - 根据上式, 我们可以对其进行理解, 通过权重将重要的位置凸显出来. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "interpreter": {
   "hash": "d59b92a5f9bdcea563f9a2c491f707fbdea763ce0564d33fd44b6f2c9829ca9d"
  },
  "kernelspec": {
   "display_name": "Python 3.7.0 64-bit ('base': conda)",
   "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.0"
  },
  "orig_nbformat": 4
 },
 "nbformat": 4,
 "nbformat_minor": 2
}