R——概率统计与模拟（六）重要性采样

作者: 生信了 | 来源:发表于2020-02-19 16:16 被阅读0次

R——概率统计与模拟（六）重要性采样
R——概率统计与模拟（二）
R——概率统计与模拟（一）
6.5.1蒙塔卡罗模拟理论概述。.where、.sqrt
R——概率统计与模拟（四）拒绝抽样
蒙特卡洛随机方法/模拟（Monte Carlo method）
数据科学家应当了解的五个统计基本概念
R——概率统计与模拟（五）彩票连号、归纳法以及二项分布
蒙特卡洛积分与重要性采样
采样的作用

原创：hxj7

本文介绍了重要性采样（Importance Sampling）并给出了示例。

本文篇幅较长，分为以下几个部分：

重要性采样是什么
重要性采样的应用示例
不同的分布 $Q$ 对结果有影响吗？
小结和附录代码

Part1：重要性采样是什么

前文《R-概率统计与模拟（三）变换均匀分布对特定分布进行抽样》和《R-概率统计与模拟（四）拒绝抽样》分别介绍了两种方法，可以根据已知的p.d.f.进行采样（抽样），使得采样得到的点符合目标分布。“重要性抽样”（Importance Sampling）不同于上述两种方法，其主要用途是通过采样计算数学期望。 而这与积分的计算有关。

通过大量随机采样计算积分

在建模过程中，我们经常要计算某个函数的期望：假设变量 $X$ 符合某个分布 $P$ ，记为 $X \sim P$ ，其 $\text{p.d.f.}$ 记为 $p(x)$ 。我们想要计算某个函数 $f$ 在分布 $P$ 下的期望 $E[f(X)]$ 。根据期望的数学定义，我们知道：
$E_{X \sim P}[f(X)] = \int_x f(x)p(x) \, \text{d}x$

在实际操作中，往往上面积分的解析解很难计算，所以可以通过大量采样的方法进行近似计算：假设根据分布 $P$ 进行大量随机采样，样本量记为 $N$ ，样本点记为 $\{\hat{x}_1, \hat{x}_2, \ldots, \hat{x}_N \}$ ，那么当 $N$ 足够大时，
$E_{X \sim P}[f(X)] = \int_x f(x)p(x) \, \text{d}x \approx \frac{1}{N} \sum_{i=1}^{N} f(\hat{x}_i)$

所以关键在于如何对分布 $P$ 进行采样。当然我们可以根据前文介绍的“变换均匀分布进行抽样”或者“拒绝抽样”的方法进行采样。但是，如前文所说，这两种方法都有局限。当遇到一个复杂的 $\text{p.d.f.}$ 时，两种方法可能都不适用。

引入另一个分布 $Q$

当无法直接对分布 $P$ 进行采样从而难以近似计算上述积分时，我们可以对上面的积分做一个变换（即“重要性采样”的关键部分），让计算变得简单可行：

我们选取另一个容易进行抽样的分布 $Q$ ，其 $\text{p.d.f.}$ 记为 $q(x)$ 。上述积分变换为：
$E_{X \sim P}[f(X)] = \int_x f(x)p(x) \, \text{d}x = \int_x f(x) \frac{p(x)}{q(x)} q(x) \, \text{dx}$
然后我们根据分布 $Q$ 进行大量随机采样，样本量记为 $N$ ，样本点记为 $\{x_1, x_2, \ldots, x_N \}$ ，那么当 $N$ 足够大时，
$E_{X \sim P}[f(X)] = \int_x f(x)p(x) \, \text{d}x = \int_x f(x) \frac{p(x)}{q(x)} q(x) \, \text{dx} \approx \frac{1}{N} \sum_{i=1}^{N} \frac{p(x_i)}{q(x_i)} f(x_i)$

相比较于原来的计算方式，变换后的式子中 $\frac{p(x_i)}{q(x_i)}$ 就是加之于 $f(x_i)$ 的权重（重要性）。

简单小结

小结一下，要计算分布 $P$ 下的数学期望 $E[f(X)]$ ，重要性采样主要分为以下两步：

选取另一个容易直接采样的分布 $Q$ ，并据此分布大量随机采样，得到一系列采样点 $\{x_1, x_2, \ldots, x_N \}$
根据上面的采样点，计算 $\frac{1}{N} \sum_{i=1}^{N} \frac{p(x_i)}{q(x_i)} f(x_i)$ 的值作为近似结果。

Part2：重要性采样的应用示例

为了验证重要性采样的可行性，本文列举了两个示例，分别是计算均匀分布和非对称拉普拉斯分布的数学期望。（注意，这两个示例本身都可以直接计算出数学期望的理论值，用它们做例子只是为了验证重要性采样是否可以得到正确的结果。）

首先看均匀分布的例子：假设变量 $X \sim \text{U}(0,1)$ ，试计算 $E[X]$ 和 $E[X^2]$ 。

很明显，上述两个期望的理论值是：
$E_T[X]=0.5, \quad E_T[X^2]=\frac{1}{3}$

其中 $T$ 表示理论值（Theoretical Value）。

接下来我们需要看看利用重要性采样的方法是否可以得到正确的结果：分布 $P$ 是简单的均匀分布，我们选取另一个简单的分布 $Q$ 来进行采样，二者的 $\text{p.d.f.}$ 如下：
$\begin{aligned} & p(x) = 1, \quad \ \ x \in [0,1] \\ & q(x) = 2x, \quad x \in [0,1] \end{aligned}$

image
我们可以利用“变换均匀分布”的方法对分布进行采样，当采样点的数量达到10万个时，结果是：

与真实值非常接近。说明采用分布 $Q$ 进行重要性采样是有效的。（这一示例的R代码见文末“代码段C1”。）

其次看一个非对称拉普拉斯分布的例子：假设变量 $X \sim \text{AL}(\mu,\sigma,p)$ ，试计算 $E[X]$ 和 $E[X^2]$ 。

本文所使用的非对称拉普拉斯分布（Asymmetric Laplace Distribution）具有如下 $\text{p.d.f.}$ ：
$p(x)=f(x;\mu,\sigma,p)=\frac{p(1-p)}{\sigma} \exp \left(-\frac{x-\mu}{\sigma} [p - I(x \le \mu)] \right)$
其中， $0<p<1$ ， $\sigma > 0$ ， $-\infty < \mu < \infty$ ， $I(\cdot)$ 是指示函数。和上面均匀分布相比，这个分布复杂了很多。根据文献[1]（Keming Yu & Jin Zhang, 2006），若 $X \sim \text{AL}(\mu,\sigma,p)$ ，则上述期望的理论值是：
$E_T[X]=\mu + \frac{\sigma(1-2p)}{p(1-p)}, \quad E_T[X^2]=\frac{\sigma^2(1-2p+2p^2)}{(1-p)^2p^2} + E_T^2[X]$

其中 $T$ 表示理论值。

我们以 $\mu=1, \sigma=0.5, p=0.2$ ，即 $X \sim \text{AL}(1, 0.5, 0.2)$ 为例，根据上面的公式，期望的理论值是：
$E_T[X]=2.875, \quad E_T[X^2]=14.90625$

接下来我们选择正态分布 $\text{N}(E_T[X], \sqrt{\text{Var}(X)})$ 作为分布 $Q$ 进行重要性采样。其中， $\text{Var}(X)=\frac{\sigma^2(1-2p+2p^2)}{(1-p)^2p^2}$ 。我们将两个分布的 $\text{p.d.f.}$ 曲线放在一起比较一下：

image
为什么选择正态分布作为分布，一是正态分布采样很方便，二是正态分布和非对称拉普拉斯分布的曲线比较相似。当采样点数量达到10万个是，结果是：

与真实值比较接近，但是没有均匀分布那个例子中的效果好（这一示例的R代码见文末“代码段C2”）。那是否换一个正态分布效果会好一点呢？这就引入一个新的问题，就是不同的分布 $Q$ 对重要性采样的结果有影响吗？

Part3：不同的分布 $Q$ 对结果有影响吗？

还以上面提到的两个例子来说明。

首先我们对 $X \sim \text{U}(0,1)$ 选用不同的分布 $Q$ ，这些分布 $Q$ 的 $\text{p.d.f.}$ 都具有如下形式：

$q(x)=ax+1-0.5a, \quad x \in [0,1]$

image
当时，就是上文中所使用的。我们对取不同的值，就会得到不同的分布。这些不同的分布的重要性采样结果如下：

image

image
图中后缀表示理论值，后缀表示分布重要性采样的结果。可以看出，不同的分布的结果都很接近于理论值，并且彼此差异不大。

我们再看看 $X \sim \text{AL}(1,0.5,0.2)$ 选用不同的分布 $Q$ 的结果

我们选用不同的正态分布作为分布 $Q$ ，这些正态分布的方差都一样，而均值以一定的位移（offset）分布在 $E_T[X]$ 的左右两侧：

image

图例中“q(x)”后面的数字是正态分布的均值相对于 $E_T[X]=2.875$ 的位移。我们取不同的位移，就会得到不同的分布。这些不同的分布 $Q$ 的重要性采样结果如下：

image

image
图中"theo"表示理论值。可以看出，不同的分布的结果是有差异的，当位移是“+5”时最接近理论值。至于为什么一些分布的结果优于另一些分布，需要后续的学习研究。
（最后提一下，当然也可以通过调整方差来得到不同的正态分布，读者可自行验证。）

Part4：小结和附录代码

“重要性抽样”（Importance Sampling）的主要用途是通过采样计算数学期望，它不是一种用于直接采样的方法。我们列举了两个示例：均匀分布 $\text{U}(0,1)$ 和非对称拉普拉斯分布 $\text{AL}(\mu,\sigma,p)$ ，通过选用对应的分布 $Q$ 采样后计算得到 $E[X]$ 和 $E[X^2]$ 的近似结果，验证了重要性采样的可行性。最后，我们讨论了采用不同的分布 $Q$ 可能对近似计算结果有影响，至于影响的原因（如何得到最优效果）有待后续学习研究。

参考

[1] "A Three-Parameter Asymmetric Laplace Distribution and Its Extension", Keming Yu & Jin Zhang, 2006. Link https://doi.org/10.1080/03610920500199018
[2] "采样方法（一）", CSDN. Link https://blog.csdn.net/Dark_Scope/article/details/70992266

最后给出文中用到的代码：

代码段C1

Unif <- function(a) {
  b <- 1 - a / 2
  Px <- function(x) 1
  Qx <- function(x) a * x + b
  ntry <- 100000
  
  set.seed(123)
  syQ0 <- runif(ntry, 0, 1)
  syQ <- (sqrt(b ^ 2 + 2 * a * syQ0) - b) / a  # sampling by q(x)
  
  onetry <- function(x) Px(x) / Qx(x) * x
  simuX <- sapply(syQ, onetry)
  EX_Q <- mean(simuX)   # simulative E(X)
  
  onetry2 <- function(x) Px(x) / Qx(x) * x ^ 2
  simuX2 <- sapply(syQ, onetry2)
  EX2_Q <- mean(simuX2)   # simulative E(X^2)
  
  c(a=a, EX_Q=EX_Q, EX2_Q=EX2_Q)
}

format_formula <- function(a) {
  b <- 1 - a / 2  # b >= 0 when a <= 2
  prefix <- suffix <- "" 
  if (a == 0) {
    return(paste0("q(x)=", b))
  } else if (a == 1) {
    if (b == 0) return("q(x)=x")
    return(paste0("q(x)=x+", b)) 
  } else {
    if (b == 0) return(paste0("q(x)=", a, "x"))
    return(paste0("q(x)=", a, "x+", b)) 
  }
}

# test
test_a <- 2   # a is the slope and a <= 2.
Unif(test_a)
tmpx <- seq(0, 1, 0.001)
tmpy <- rep(1, length(tmpx))
npoint <- 1000
plot(tmpx, tmpy, col=1, xlab="x", ylab="y", main="p.d.f. curves", cex=.1, 
     ylim=c(0, test_a / 2 + 1.1))
curve(test_a * x + (1 - test_a / 2), 0, 1, n=npoint, add=T, col=2)
legend("topleft", legend=c("p(x)=1", format_formula(test_a)), 
       col=1:2, lty=1, bty="n", cex=.8, ncol=2)

# different values of a.
all_a <- c(0.1, 0.2, 0.5, 1, 2)  # a is the slope and a <= 2.
res <- as.data.frame(t(sapply(all_a, Unif)))

# plot result of different values of a
res$EX_T <- 0.5   # theoretical E(X)
res$EX2_T <- 1 / 3  # theoretical E(X^2)
res$a <- factor(res$a, levels=res$a)
res
library(dplyr)
library(tidyr)
library(ggplot2)
res %>%
  gather(Expectation, Value, -a) %>%
  ggplot(aes(x=a, y=Value, color=Expectation, group=Expectation)) +
  geom_point() +
  geom_line() +
  labs(x="Slope")

# plot p.d.f. curves
tmpx <- seq(0, 1, 0.001)
tmpy <- rep(1, length(tmpx))
npoint <- 1000
plot(tmpx, tmpy, col=1, xlab="x", ylab="y", main="p.d.f. curves", cex=.1, 
     ylim=c(0, max(all_a) / 2 + 1.1))
i <- 2
for (a in all_a) {
  curve(a * x + (1 - a / 2), 0, 1, n=npoint, add=T, col=i)
  i <- i + 1
}
legend("topleft", legend=c("p(x)", sapply(all_a, format_formula)),
       col=1:(length(all_a) + 1), lty=1, bty="n", cex=.8, ncol=2)

代码段C2

# p.d.f. of the asymmetric Laplace distribution.
u <- 1
sg <- 0.5
p <- 0.2
ALx <- function(x, u, sg, p) {
  if (x >= u) {
    return(p * (1 - p) / sg * exp(-p * (x - u) / sg))
  } else {
    return(p * (1 - p) / sg * exp((1 - p) * (x - u) / sg))
  }
}
Px <- function(x) ALx(x, u=u, sg=sg, p=p)

# theoretical value
EX_T <- C1 <- sg * (1 - 2 * p) / (p * (1 - p)) + u   # theoretical E(X)
C2 <- sg / (p * (1 - p)) * sqrt((1 - 2 * p + 2 * p ^ 2))
EX2_T <- C2 ^ 2 + C1 ^ 2  # theoretical E(X^2)

# sampling/simulation parameters.
ntry <- 100000

# importance sampling
IS <- function(offset) {
  Qx <- function(x) dnorm(x, C1 + offset, C2)  # use normal distribution as q(x)
  set.seed(123)
  sy <- rnorm(ntry, C1 + offset, C2) # sampling by q(x)
  
  onetryQ <- function(x) Px(x) / Qx(x) * x
  simu_Q <- sapply(sy, onetryQ)
  EX_Q <- mean(simu_Q)  # simulative E(X) by importance sampling
  
  onetryQ2 <- function(x) Px(x) / Qx(x) * x ^ 2
  simu_Q2 <- sapply(sy, onetryQ2)
  EX2_Q <- mean(simu_Q2) # simulative E(X^2) by importance sampling
  
  c(Group=offset, EX=EX_Q, EX2=EX2_Q)
}

# test
test_offset <- 0   # offset from theoretical E(X)
IS(test_offset)
npoints <- 10000
xbegin <- -10
xend <- 15
plotPx <- function(x) sapply(x, Px)
curve(plotPx, xbegin, xend, n=npoints, col=1, xlab="x", ylab="y", main="p.d.f. curves")
Qx <- function(x) dnorm(x, C1, C2)
plotQx <- function(x) sapply(x, Qx)
curve(plotQx, xbegin, xend, n=npoints, add=T, col=2)
legend("topleft", legend=c("p(x)", "q(x)"), 
       col=1:2, lty=1, cex=.8, bty="n", ncol=2)

# different values of offset
all_offset <- c(-5, -2, 0, 2, 5, 10)
offset_labels <- c("-5", "-2", "0", "+2", "+5", "+10")
res <- as.data.frame(t(sapply(all_offset, IS)))
res$Group <- offset_labels
res <- rbind(res, c("theo", EX_T, EX2_T))
res$EX <- as.numeric(res$EX)
res$EX2 <- as.numeric(res$EX2)
res$Group <- factor(res$Group, levels=c(offset_labels, "theo"))
res

# plot the result of different values of offset.
library(dplyr)
library(tidyr)
library(ggplot2)
res %>%
  gather(Expectation, Value, -Group) %>%
  ggplot(aes(x=Expectation, y=Value, fill=Group)) +
  geom_bar(stat="identity", position="dodge")

# plot the asymmetric Laplace distribution.
npoints <- 10000
xbegin <- -10
xend <- 15
plotPx <- function(x) sapply(x, Px)
curve(plotPx, xbegin, xend, n=npoints, col=1, xlab="x", ylab="y", main="p.d.f. curves")
i <- 2
for (offset in all_offset) {
  Qx <- function(x) dnorm(x, C1 + offset, C2)
  plotQx <- function(x) sapply(x, Qx)
  curve(plotQx, xbegin, xend, n=npoints, add=T, col=i)
  i <- i + 1
}
legend("topleft", legend=c("p(x)", paste0("q(x),", offset_labels)), 
       col=1:(i - 1), lty=1, cex=.8, bty="n", ncol=2, x.intersp=.5, text.width=2.5)

（公众号：生信了）

R——概率统计与模拟（六）重要性采样
原创：hxj7 本文介绍了重要性采样（Importance Sampling）并给出了示例。本文篇幅较长，分为以...
R——概率统计与模拟（二）
原创：hxj7 本文继续介绍一些和概率统计相关的模拟。前文《R-概率统计与模拟》介绍了一些用 R 进行概率模拟的...
R——概率统计与模拟（一）
原创：hxj7 本文记录了三个概率统计相关的小题目，以回顾一些概率统计的知识。正如笔者在前文《公众号一岁啦》中所...
6.5.1蒙塔卡罗模拟理论概述。.where、.sqrt
蒙塔卡罗模拟：通过随机采样来找到问题的近似答案 '蒙特卡罗模拟方法，又称，‘随机抽样’或‘统计实验方法，是以概率和...
R——概率统计与模拟（四）拒绝抽样
原创：hxj7 本文介绍了拒绝抽样（Reject Sampling）。前文《R-概率统计与模拟（三）变换均匀分布...
蒙特卡洛随机方法/模拟（Monte Carlo method）
蒙特卡洛随机方法，即统计模拟方法，是一类以概率统计理论为指导的数值计算方法。本质上是用部分估计整体，采样越多，则越...
数据科学家应当了解的五个统计基本概念
摘要：本文讲述了数据科学家应当了解的五个统计基本概念：统计特征、概率分布、降维、过采样/欠采样、贝叶斯统计从高...
R——概率统计与模拟（五）彩票连号、归纳法以及二项分布
原创：hxj7 《R-概率统计与模拟》的系列文章已经写了四篇。其中有经过整理，内容相对全面，并且有内在联系的文章；...
蒙特卡洛积分与重要性采样
重要性采样在强化学习有着重要作用,它是蒙特卡洛积分的一种采样策略. 目录概率论基础蒙特卡洛积分重要性采样参...
采样的作用
定义：采样本质上是对随机现象的模拟，根据给定的概率分布，来模拟产生一个对应的随机事件。采样可以让人们对随机事件及其...

R——概率统计与模拟（六）重要性采样

本文篇幅较长，分为以下几个部分：

Part1：重要性采样是什么

通过大量随机采样计算积分

引入另一个分布 $Q$

简单小结

Part2：重要性采样的应用示例

首先看均匀分布的例子：假设变量 $X \sim \text{U}(0,1)$ ，试计算 $E[X]$ 和 $E[X^2]$ 。

其次看一个非对称拉普拉斯分布的例子：假设变量 $X \sim \text{AL}(\mu,\sigma,p)$ ，试计算 $E[X]$ 和 $E[X^2]$ 。

Part3：不同的分布 $Q$ 对结果有影响吗？

首先我们对 $X \sim \text{U}(0,1)$ 选用不同的分布 $Q$ ，这些分布 $Q$ 的 $\text{p.d.f.}$ 都具有如下形式：

我们再看看 $X \sim \text{AL}(1,0.5,0.2)$ 选用不同的分布 $Q$ 的结果

Part4：小结和附录代码

参考

代码段C1

代码段C2

相关文章