重要性采样(Importance Sampling)是一种用于估计期望值的统计方法,特别适用于从难以直接采样的分布中获取样本的情况。其核心思想是通过一个易于采样的提议分布(proposal distribution)来近似目标分布(target distribution),从而计算期望值。
重要性采样的本质是计算期望值,在《概率论》这门学科中,计算期望值是最核心的基础操作,重要性采样从另外一个角度告诉人们如何计算期望值的方式,它通过从一个易于采样的分布中抽取样本,并对这些样本进行加权,来估计目标分布的期望值。
目标分布:我们想要从中采样的分布 $p(x)$,通常是难以直接采样的。
提议分布:一个易于采样的分布$q(x)$,用于近似目标分布。
重要性权重:用于调整从提议分布中采样的样本,使其能够反映目标分布的特性。权重定义为$w(x_i) = \frac{p(x_i)}{q(x_i)}$。
假设我们想要估计函数 $f(x)$ 在目标分布 $p(x)$ 下的期望值:
$$\mathbb{E}_{p(x)}[f(x)] = \int f(x) p(x) \, dx$$
通过重要性采样,可以将其转化为在提议分布$q(x)$下的期望值:
$$\mathbb{E}_{p(x)}[f(x)] = \int f(x) \frac{p(x)}{q(x)} q(x) \,dx$$
$$= \mathbb{E}_{q(x)}[f(x) \frac{p(x)}{q(x)}]$$
因此,我们可以从 $q(x)$ 中采样,并计算加权平均: $$\hat{\mu} = \frac{1}{N} \sum_{i=1}^N f(x_i) w(x_i)$$
其中,$w(x_i) = \frac{p(x_i)}{q(x_i)}$
能够处理难以直接采样的目标分布。
通过选择合适的提议分布,可以提高估计的效率和准确性。
如果提议分布 $q(x)$ 与目标分布$p(x)$ 差异较大,可能导致权重分布不均,估计方差增大。
需要计算重要性权重,可能增加计算复杂度。
重要性采样广泛应用于蒙特卡洛方法、贝叶斯推断、强化学习等领域,尤其是在高维空间或复杂分布中。
既然目标分布难以采样,怎么在计算权重的时候又出现了目标分布呢?
这是一个非常好的问题!重要性采样中的核心矛盾在于:目标分布 $p(x)$ 难以直接采样,但我们仍然需要知道它的形式(或至少知道它的未归一化形式)来计算权重。
重要性采样是一种通过改变采样分布来提高估计效率的方法。以下是一个通俗易懂的例子(估计一个稀有事件的概率):
假设你有一个不均匀的骰子,掷出6的概率非常低(比如0.1%),你想估计这个概率。如果直接用蒙特卡洛方法,可能需要掷很多次骰子才能观察到几次6,效率很低。
在此情况下,可以应用重要性采样,总共分为下面三个过程:
1、改变采样分布:你可以设计一个“偏向”骰子,使得掷出6的概率更高(比如10%)。这样,你更容易观察到6。
2、加权调整:每次掷出6时,记录结果并乘以一个权重,这个权重是原始分布与改变后分布的概率比。例如,如果原始概率是0.1%,改变后是10%,那么权重就是0.1% / 10% = 0.01。
3、计算估计值:通过多次掷骰,记录加权结果并取平均,得到更准确的估计。
总之,重要性采样通过增加稀有事件的采样频率,并用权重调整结果,从而在较少样本下获得更精确的估计。这个方法在金融风险评估、物理模拟等领域有广泛应用。
下面是使用 PyTorch 实现重要性采样的代码。我们将使用 PyTorch 的张量操作和随机数生成功能来实现上述逻辑。
import torch
# 目标分布 p(x) = exp(-x^2 / 2) / sqrt(2 * pi)
def target_distribution(x):
return torch.exp(-x**2 / 2) / torch.sqrt(2 * torch.tensor(torch.pi))
# 提议分布 q(x) = 1 / (1 + x^2) / pi
def proposal_distribution(x):
return 1 / (1 + x**2) / torch.tensor(torch.pi)
# 从提议分布中采样(柯西分布)
def sample_from_proposal(n_samples):
return torch.tan(torch.tensor(torch.pi) * (torch.rand(n_samples) - 0.5))
# 重要性采样
def importance_sampling(n_samples):
samples = sample_from_proposal(n_samples)
weights = target_distribution(samples) / proposal_distribution(samples)
return samples, weights
# 估计期望值
def estimate_expectation(n_samples):
samples, weights = importance_sampling(n_samples)
expectation = torch.mean(samples * weights)
return expectation
# 参数设置
n_samples = 10000
# 估计期望值
estimated_expectation = estimate_expectation(n_samples)
print(f"Estimated Expectation: {estimated_expectation.item()}")
target_distribution(x) 是标准正态分布的概率密度函数。
使用 PyTorch 的 torch.exp 和 torch.sqrt 实现。
proposal_distribution(x) 是柯西分布的概率密度函数。
柯西分布的公式为 $q(x)= \frac{ π(1+x^2)}{1}$
柯西分布可以通过均匀分布生成:$x = \tan(\pi (u - 0.5))$,其中$u \sim \text{Uniform}(0,1)$。
使用 torch.rand 生成均匀分布,并通过变换得到柯西分布的样本。
计算每个样本的权重:$w(x_i) = \frac{p(x_i)}{q(x_i)}$
使用 PyTorch 的张量操作计算权重。
使用加权样本计算期望值:$\mathbb{E}[x] \sim \tfrac{1}{N} \sum_{i=1}^N x_i w(x_i)$
使用 torch.mean 计算加权平均值。
运行代码后,你会得到一个估计的期望值。由于目标分布是标准正态分布,期望值应该接近 0。例如:
Estimated Expectation: 0.0123
如果目标分布完全未知,重要性采样的直接应用会变得困难,因为我们需要知道目标分布的概率密度函数(PDF)来计算权重。然而,即使目标分布未知,仍然有一些方法可以尝试解决这个问题。以下是几种可能的解决方案:
如果目标分布未知,但我们可以从目标分布中采样(例如通过实验或模拟),则可以通过以下步骤近似目标分布:
步骤 1:从目标分布中收集一些样本。
步骤 2:使用这些样本来拟合一个近似的概率分布(例如高斯混合模型、核密度估计等)。
步骤 3:将拟合的分布作为目标分布的近似,然后进行重要性采样。
代码示例:使用核密度估计(KDE)近似目标分布
import torch
from sklearn.neighbors import KernelDensity
# 假设我们有一些从目标分布中采样的数据
target_samples = torch.randn(1000) # 例如,目标分布是标准正态分布
# 使用核密度估计(KDE)拟合目标分布
kde = KernelDensity(kernel="gaussian", bandwidth=0.5).fit(target_samples.reshape(-1, 1))
# 定义近似的目标分布
def approximate_target_distribution(x):
log_prob = kde.score_samples(x.reshape(-1, 1))
return torch.tensor(np.exp(log_prob))
# 提议分布(例如柯西分布)
def proposal_distribution(x):
return 1 / (1 + x**2) / torch.tensor(torch.pi)
# 从提议分布中采样
def sample_from_proposal(n_samples):
return torch.tan(torch.tensor(torch.pi) * (torch.rand(n_samples) - 0.5))
# 重要性采样
def importance_sampling(n_samples):
samples = sample_from_proposal(n_samples)
weights = approximate_target_distribution(samples) / proposal_distribution(samples)
return samples, weights
# 估计期望值
def estimate_expectation(n_samples):
samples, weights = importance_sampling(n_samples)
expectation = torch.mean(samples * weights)
return expectation
# 参数设置
n_samples = 10000
# 估计期望值
estimated_expectation = estimate_expectation(n_samples)
print(f"Estimated Expectation: {estimated_expectation.item()}")
如果目标分布完全未知且无法采样,但仍然可以评估某个函数在目标分布下的期望值(例如通过实验或模拟),则可以使用无模型方法,例如:
蒙特卡洛方法:直接通过实验或模拟生成样本,计算样本的平均值。
强化学习中的离策略方法:在强化学习中,即使目标策略未知,也可以通过行为策略采样并使用重要性采样来估计目标策略的价值函数。
自适应重要性采样是一种迭代方法,通过逐步改进提议分布来逼近目标分布。即使目标分布未知,也可以通过以下步骤实现:
步骤 1:初始化一个提议分布 q(x)。
步骤 2:从提议分布中采样,并根据样本调整提议分布(例如通过最大化似然或最小化方差)。
步骤 3:重复步骤 2,直到提议分布足够接近目标分布。
代码示例:自适应重要性采样
import torch
import numpy as np
# 假设目标分布是某种复杂分布,我们无法直接知道它的形式
# 但我们有一个黑箱函数可以计算目标分布的概率密度
def target_distribution(x):
return torch.exp(-x**2 / 2) * (1 + torch.sin(x * 2)) # 示例目标分布
# 初始提议分布(例如高斯分布)
def proposal_distribution(x, mu, sigma):
return torch.exp(-0.5 * ((x - mu) / sigma)**2) / (sigma * torch.sqrt(2 * torch.tensor(torch.pi)))
# 从提议分布中采样
def sample_from_proposal(n_samples, mu, sigma):
return torch.normal(mu, sigma, size=(n_samples,))
# 自适应重要性采样
def adaptive_importance_sampling(n_samples, n_iterations):
mu = 0.0 # 初始均值
sigma = 1.0 # 初始标准差
for _ in range(n_iterations):
samples = sample_from_proposal(n_samples, mu, sigma)
weights = target_distribution(samples) / proposal_distribution(samples, mu, sigma)
# 更新提议分布的参数
mu = torch.sum(samples * weights) / torch.sum(weights)
sigma = torch.sqrt(torch.sum((samples - mu)**2 * weights) / torch.sum(weights))
return samples, weights
# 估计期望值
def estimate_expectation(n_samples, n_iterations):
samples, weights = adaptive_importance_sampling(n_samples, n_iterations)
expectation = torch.mean(samples * weights)
return expectation
# 参数设置
n_samples = 10000
n_iterations = 10
# 估计期望值
estimated_expectation = estimate_expectation(n_samples, n_iterations)
print(f"Estimated Expectation: {estimated_expectation.item()}")
如果目标分布完全未知,但可以评估某个函数在目标分布下的期望值(例如通过实验或模拟),则可以使用黑箱优化方法(如贝叶斯优化)来直接优化目标函数。
当目标分布完全未知时,可以通过以下方法解决重要性采样问题:
使用样本数据拟合目标分布的近似分布。
使用无模型方法直接估计期望值。
使用自适应重要性采样逐步逼近目标分布。
使用黑箱优化方法直接优化目标函数。
选择哪种方法取决于具体问题的性质以及可用的信息和资源。