重要性采样求泊松分布期望

我也不是数值模拟或者统计物理的学生。所以蒙特卡洛之类的概念性的东西也说的不好。下面是我的统计物理的期末作业。关于用重要性分布求离散函数泊松分布的数学期望的。
我觉得这个挺有意思的，所以也花了点时间加上chatgpt折腾了一下，提供一个参考吧，如果有不对的地方，欢迎讨论
具体直接看资源绑定的pdf吧。
或者下面分块放到jupyter notebook里面会更加清楚一点

蒙特卡罗采样

使用Metropolis 算法这样我们可以对任意的采样函数进行采样

import numpy as np
import matplotlib.pyplot as plt

目标分布 g(x) (未归一化)

def g(x):
return x**3

Metropolis 算法 (仅生成整数样本)

def metropolis_sampling(g, x0, delta, num_samples, a, b):
samples = []
x = x0
for _ in range(num_samples):
y = x + np.random.randint(-delta, delta + 1) # 生成整数的提议点
if y < a or y > b:
samples.append(x)
else:
acceptance_ratio = g(y) / g(x)
if np.random.rand() < acceptance_ratio:
x = y
samples.append(x)
return np.array(samples)

设置参数

x0 = 5 # 初始点
delta = 1 # 调整参数
num_samples = 100000 # 样本数
a = 0 # 区间下限
b = 10 # 区间上限

运行 Metropolis 算法

samples = metropolis_sampling(g, x0, delta, num_samples, a, b)

绘制样本分布

plt.hist(samples, bins=np.arange(a, b+1)-0.5, density=True, alpha=0.75, edgecolor=‘black’, label=‘Sampled Distribution’)

绘制目标分布（归一化）

x = np.arange(a, b+1)
y = g(x)
y = y / np.sum(y) # 归一化
plt.plot(x, y, ‘r-’, lw=2, label=‘Target Distribution $g (x)$ ’)

plt.xlabel(‘x’)
plt.ylabel(‘Density’)
plt.title(‘Metropolis Sampling (Integer Samples)’)
plt.legend()
plt.grid(True)
plt.show()

可以看到Metropolis 算法很好的对x^2这个函数进行了采样。因为泊松方程只能取整数，然后我才对蒙特卡洛采样程序进行修改，让他样本的结果只为整数

使用正态分布作为重要性采样

重要性采样函数力求与目标函数接近，正态函数与泊松函数的接近程度更高一点

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import norm
from scipy.stats import poisson

目标分布 g(x) (未归一化的正态分布)

def g(x, mu=5, sigma=1):
return np.exp(-0.5 * ((x - mu) / sigma) ** 2)

Metropolis 算法 (仅生成整数样本)

def metropolis_sampling(g, x0, delta, num_samples, a, b, mu=5, sigma=1):
samples = []
x = x0
for _ in range(num_samples):
y = x + np.random.randint(-delta, delta + 1) # 生成整数的提议点
if y < a or y > b:
samples.append(x)
else:
acceptance_ratio = g(y, mu, sigma) / g(x, mu, sigma)
if np.random.rand() < acceptance_ratio:
x = y
samples.append(x)
return np.array(samples)

设置参数

x0 = 5 # 初始点
delta = 1 # 调整参数
num_samples = 100000 # 样本数
a = 0 # 区间下限
b = 16 # 区间上限
mu = 5 # 正态分布的均值
sigma = 1 # 正态分布的标准差

运行 Metropolis 算法

samples = metropolis_sampling(g, x0, delta, num_samples, a, b, mu, sigma)

绘制样本分布

plt.hist(samples, bins=np.arange(a, b+1)-0.5, density=True, alpha=0.75, edgecolor=‘black’, label=‘Sampled Distribution’)

绘制目标分布（归一化）

x = np.arange(a, b+1)
y = g(x, mu, sigma)
y = y / np.sum(y) # 归一化
plt.plot(x, y, ‘r-’, lw=2, label=‘Target Distribution $g (x)$ ’)

plt.xlabel(‘x’)
plt.ylabel(‘Density’)
plt.title(‘Metropolis Sampling (Integer Samples, Normal Distribution)’)
plt.legend()
plt.grid(True)
plt.show()
#对正态样本100000的统计图，采样的样本趋于原来的正态函数

计算权重 w(x) = P(x) / q(x)

lambda_param = 8
weights = poisson.pmf(samples, lambda_param) / norm.pdf(samples, mu, sigma)

根据权重计算泊松分布的期望值

expectation = np.mean(weights * samples)
print(‘期望值%f’ %expectation)
plt.plot(samples,norm.pdf(samples, mu, sigma),‘*’,label=‘norm’)
plt.plot(samples,poisson.pmf(samples, lambda_param),‘.’,label=‘poisson’)
plt.legend()
plt.show()

这个泊松分布的结果只能用糟糕了形容。从上图采样函数norm和目标函数poisson也可以看出来，这个正态函数的参数取的并不好。
修改一下采样函数的均值另外对于正态分布的标准差也进行一定的修改

设置参数

x0 = 8 # 初始点
delta = 1 # 调整参数
num_samples = 100000 # 样本数
a = 0 # 区间下限
b = 16 # 区间上限
mu = 8 # 正态分布的均值

运行 Metropolis 算法

for sigma in range(1,5):
samples = metropolis_sampling(g, x0, delta, num_samples, a, b, mu, sigma)
# 计算权重 w(x) = P(x) / q(x)
lambda_param = 8
weights = poisson.pmf(samples, lambda_param) / norm.pdf(samples, mu, sigma)

# 根据权重计算泊松分布的期望值
expectation = np.mean(weights * samples)print('sigma=%d' %sigma)
print(f'泊松分布的期望值: {expectation:.5f}')
plt.plot(samples,norm.pdf(samples, mu, sigma),'*',label='norm')   
plt.plot(samples,poisson.pmf(samples, lambda_param),'.',label='poisson')
plt.legend()
plt.show()
print('\n')

至此发现当正态函数的sigma取到3，mu取8的正态函数分布采样能够更加准确的表示泊松函数的形态

设置参数

from IPython.display import display, Math
x0 = 8 # 初始点
delta = 1 # 调整参数
a = 0 # 区间下限
b = 16 # 区间上限
mu = 8 # 正态分布的均值
sigma=3
print(‘误差值计算公式如下’)
formula = r"\frac{|I - I^{{\text{true}}|}{I}{\text{true}}}"
display(Math(formula))
print(‘\n’)
large_list = [10**i for i in range(2, 10)] #生成不同的样本数

for num_samples in large_list:
samples = metropolis_sampling(g, x0, delta, num_samples, a, b, mu, sigma)
# 计算权重 w(x) = P(x) / q(x)
lambda_param = 8
weights = poisson.pmf(samples, lambda_param) / norm.pdf(samples, mu, sigma)

# 根据权重计算泊松分布的期望值
expectation = np.mean(weights * samples)print('样本数=%d' %num_samples)
print(f'泊松分布的期望值: {expectation:.5f}')
print('误差值%s:' %(abs(round(expectation, 5)-lambda_param)/lambda_param))plt.plot(samples,norm.pdf(samples, mu, sigma),'*',label='norm')   
plt.plot(samples,poisson.pmf(samples, lambda_param),'.',label='poisson')
plt.legend()
plt.show()
print('\n')