r语言网站开发,搜索引擎优化的英文,中国制造网效果怎么样,男女做爰视频网站在线我也不是数值模拟或者统计物理的学生。所以蒙特卡洛之类的概念性的东西也说的不好。下面是我的统计物理的期末作业。关于用重要性分布求离散函数泊松分布的数学期望的。 我觉得这个挺有意思的#xff0c;所以也花了点时间加上chatgpt折腾了一下#xff0c;提供一个参考吧所以也花了点时间加上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, binsnp.arange(a, b1)-0.5, densityTrue, alpha0.75, edgecolor‘black’, label‘Sampled Distribution’)
绘制目标分布归一化
x np.arange(a, b1) y g(x) y y / np.sum(y) # 归一化 plt.plot(x, y, ‘r-’, lw2, label‘Target Distribution g ( x ) g(x) 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, mu5, sigma1): return np.exp(-0.5 * ((x - mu) / sigma) ** 2)
Metropolis 算法 (仅生成整数样本)
def metropolis_sampling(g, x0, delta, num_samples, a, b, mu5, sigma1): 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, binsnp.arange(a, b1)-0.5, densityTrue, alpha0.75, edgecolor‘black’, label‘Sampled Distribution’)
绘制目标分布归一化
x np.arange(a, b1) y g(x, mu, sigma) y y / np.sum(y) # 归一化 plt.plot(x, y, ‘r-’, lw2, label‘Target Distribution g ( x ) g(x) 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),*,labelnorm)
plt.plot(samples,poisson.pmf(samples, lambda_param),.,labelpoisson)
plt.legend()
plt.show()
print(\n)至此发现当正态函数的sigma取到3mu取8的正态函数分布采样能够更加准确的表示泊松函数的形态
设置参数
from IPython.display import display, Math x0 8 # 初始点 delta 1 # 调整参数 a 0 # 区间下限 b 16 # 区间上限 mu 8 # 正态分布的均值 sigma3 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),*,labelnorm)
plt.plot(samples,poisson.pmf(samples, lambda_param),.,labelpoisson)
plt.legend()
plt.show()
print(\n)总结
1泊松分布的期望为参数λ8 2重要性采样中采样函数和目标函数越解决期望越准确 3根据大数定律样本数越多估算的期望越准确 4目标函数如果只能取整数那么样本也只能取整数
