积分倒数的无偏估计

积分倒数的无偏估计

问题引出

  • 已知 \(X\sim U(1,2)\)\(p(x)=1\)

\[ g(x)=x^2 \]

  • 求值:\(E\left[g(X)\right]\)\(\dfrac{1}{E\left[g(X)\right]}\)\(E\left[\dfrac{1}{g(X)}\right]\)

解析结果

\[ E[g(X)]=\int_{1}^{2}x^2p(x)\;\mathrm{d}x=\dfrac{x^3}{3}\Bigg\vert_{1}^{2}=\dfrac{7}{3} \]

\[ \dfrac{1}{E\left[g(X)\right]}=\dfrac{3}{7}\approx0.4285714 \]

\[ E\left[\dfrac{1}{g(X)}\right]=\int_{1}^{2}\dfrac{1}{x^2}p(x)\;\mathrm{d}x=-\dfrac{1}{x}\Bigg\vert_{1}^{2}=\dfrac{1}{2} \]

  • 这里我们可以看到 \(E\left[g(X)\right]\ne E\left[\dfrac{1}{g(X)}\right]\)

无偏估计

  • \(E\left[g(X)\right]\) 的无偏估计 \(\hat{g}_1\) 如下

\[ \hat{g}_1=\dfrac{1}{N}\sum_{i=1}^{N}\dfrac{x_i^2p(x_i)}{q(x_i)} \]

  • \(q(x)\) 可以和 \(p(x)\) 相同,但是这里取一个不同的值,方便看结果
  • 例如:\(q(x)\propto x\)

\[ q(x)=x\Big/\int_{1}^{2}xp(x)\;\mathrm{d}x=\dfrac{2}{3}x \]

  • 如下估计 \(\hat{g}_2\) 并不是 \(\dfrac{1}{E\left[g(X)\right]}\) 的无偏估计,而是 \(E\left[\dfrac{1}{g(X)}\right]\) 的无偏估计

\[ \hat{g}_2=\dfrac{1}{N}\sum_{i=1}^{N}\dfrac{1}{x_i^2}\cdot\dfrac{p(x_i)}{q(x_i)} \]

采样细节

  • 如何从均匀分布变化:\(U(1,2)\)
    • 目标 CDF:\(F(x)=\dfrac{x^2-1}{3}\)
      • \(F(x)=\text{Pr}(X\le x)=\text{Pr}(f(U)\le x)=\text{Pr}(U\le f^{-1}(x))=f^{-1}(x)-1\)
    • \(f(u(x))=\sqrt{3u(x)-2}\)

问题

  • 我们无法得到 \(E\left[g(X)\right]\) 的真值,如何从 \(E\left[g(X)\right]\) 的无偏估计获取 \(E\left[\dfrac{1}{g(X)}\right]\) 的无偏估计?

模拟

  • 若干样本点:\(x_1,\cdots,x_{N_1},x_{N_1+1},\cdots,x_{N_2},x_{N_2+1},\cdots,x_{N_M}\)
  • \(\hat{g}_1\) 不是 \(\alpha\) 的无偏估计,\(\hat{g}_2\)不是 \(\alpha\) 的无偏估计
    • \(E[\hat{g}_1]\ne\alpha,E[\hat{g}_2]\ne\alpha\)

    • 需要无偏吗?无偏是为了什么?

    • 可以这么计算,但是这样并不是无偏估计

    • \({\color{red}M\to\infty}\),此时得到的结果才是正确的

\[ \hat{g}_1=\dfrac{1}{\dfrac{1}{N_1}\sum_{i=1}^{N_1}f(x_i)/p(x_i)} \]

\[ \hat{g}_2=\dfrac{1}{M}\left(\dfrac{1}{\dfrac{1}{N_1}\sum_{i=1}^{N_1}f(x_i)/p(x_i)}+\cdots+\dfrac{1}{\dfrac{1}{N_M}\sum_{i=N+1}^{N_M}f(x_i)/p(x_i)}\right) \]

  • 数值模拟
    • \(N=1,M\to\infty\):结果并不是真实值
    • \(N\to\infty,M=1\):结果是真实值

类光线追踪方法

  • 通过泰勒展开的方式转化为 PT 的形式,加以计算
  • 是无偏估计

其他方式

  • 是无偏估计

  • 和上面其实类似

  • 从概率分布 \(p(x)\) 中随机采样 \(x_i\) ,以 \(\dfrac{f(x_i)}{p(x_i)}\) 的概率终止 ,以 \(1-\dfrac{f(x_i)}{p(x_i)}\) 的概率继续

    • 要求:\(\dfrac{f(x_i)}{p(x_i)}\in(0,1),\forall x_i\)
  • 采样次数 \(t\) 的期望 \(E[t]\) 就是 \(\dfrac{1}{E\left[g(X)\right]}\) 的值

\[ E[t]=1+E[t]\int\left(1-\dfrac{f(x)}{p(x)}\right)\cdot p(x)\;\mathrm{d}x \]

\[ E[t]=\dfrac{1}{\int f(x)\cdot p(x)\;\mathrm{d}x} \]

  • 另外的理解,构造随机变量 \(X\) 如下,求期望得到上面式子
    • \(X_i\) 为和 \(X\) 独立同分布的随机变量
    • \(Y\) 是满足分布 \(p(\cdot)\) 的随机变量

\[ X=1+\left(1-\dfrac{f(Y)}{p(Y)}\right)X_i \]

  • 可以让 \(f(x)\) 缩放到 \((0,1)\),估计 \(1\Big/\int k\cdot f(x)\cdot p(x)\;\mathrm{d}x\),然后再将结果 \(\times k\)
    • 这样方便控制 \(p(x)\)

代码

其他

  • 如何估计

\[ \int\dfrac{f(x)}{\int p(x\mid y)\;\mathrm{d}y}\;\mathrm{d}x \]

\[ \int \left(f(x)\int p(x\mid y)\;\mathrm{d}y\right)\;\mathrm{d}x \]