当前位置：首页 > news >正文

利用MCMC 获得泊松分布

news 2025/8/26 7:40:00

$P_n=\frac{\lambda^nexp(-\lambda)}{n!},n=0,1,2,...$

$X\sim P(\lambda),E(X)=\sqrt{D(X)}=\lambda$

写出概率流方程如下

        if state == 0:            if np.random.random() <= min([Lambda/2, 1]):state = 1else:passelif state == 1:if choose_prob_state[i] <= 0.5:#选择 1 -> 0，此时的接受概率为min[2/Lambda, 1]if np.random.random() <= min([2/Lambda, 1]):state = 0else:passelse:#选择 1 -> 2，此时接受概率为 min[Lambda/(n+1), 1]if np.random.random() <= min([Lambda/(state+1), 1]):state = 2else:passelif state >= 2:if choose_prob_state[i] <= 0.5:#选择 n -> n+1，此时接受概率为 min[Lambda/(n+1), 1]if np.random.random() <= min([Lambda/(state+1), 1]):state = state + 1else:passelse:#选择 n+1 > n，此时接受概率为 min[(n+1)/Lambda, 1]if np.random.random() <= min([(state)/Lambda, 1]):state = state - 1else:pass

blocking 方法

def block_averages(data, block_size):num_blocks = len(data) // block_sizeblocks = data[:num_blocks*block_size].reshape(num_blocks, block_size)block_avgs = blocks.mean(axis=1)return block_avgsblock_mean = []
block_std  = []for i in range(1, 201):block_size = 5 * iblock_avgs = block_averages(results, block_size)mean_estimate = np.mean(block_avgs)standard_error = np.std(block_avgs, ddof=1) / np.sqrt(len(block_avgs))block_mean.append(mean_estimate)block_std.append(standard_error)

Lambda = 1 生成效果

average time: 1.072e-06
ave: 0.9996688
std: 1.00027000870093
(array([3.681131e+06, 3.678446e+06, 1.837276e+06, 6.127200e+05,
1.533770e+05, 3.116400e+04, 5.095000e+03, 7.020000e+02,
8.300000e+01, 6.000000e+00]), array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10.]), <BarContainer object of 10 artists>)

blocking method

随着block 增大稳定效果显著

Lambda = 7

average time: 1.153e-06
ave: 7.0095212
std: 2.6496322285839153
(array([9.062000e+03, 6.352700e+04, 2.216480e+05, 5.190980e+05,
9.097340e+05, 1.274978e+06, 1.487161e+06, 1.487430e+06,
1.304976e+06, 1.016897e+06, 7.126600e+05, 4.541560e+05,
2.646540e+05, 1.432550e+05, 7.228000e+04, 3.374700e+04,
1.474600e+04, 6.073000e+03, 2.455000e+03, 9.640000e+02,
3.790000e+02, 9.900000e+01, 1.700000e+01, 4.000000e+00]), array([ 0., 1., 2., 3., 4., 5., 6., 7., 8., 9., 10., 11., 12.,
13., 14., 15., 16., 17., 18., 19., 20., 21., 22., 23., 24.]), <BarContainer object of 24 artists>)

完整代码

import numpy as np
import matplotlib.pyplot as plt
np.random.seed(20)
import copy
import time##pn = \lambda^n * exp(-\lambda)/n!def poidis(Lambda, num, init=0):random_list = np.zeros(num)state = initmax_state = initrandom_list[0] = statechoose_prob_state = np.random.random(num)for i in range(1, num):if state == 0:            if np.random.random() <= min([Lambda/2, 1]):state = 1else:passelif state == 1:if choose_prob_state[i] <= 0.5:#选择 1 -> 0，此时的接受概率为min[2/Lambda, 1]if np.random.random() <= min([2/Lambda, 1]):state = 0else:passelse:#选择 1 -> 2，此时接受概率为 min[Lambda/(n+1), 1]if np.random.random() <= min([Lambda/(state+1), 1]):state = 2else:passelif state >= 2:if choose_prob_state[i] <= 0.5:#选择 n -> n+1，此时接受概率为 min[Lambda/(n+1), 1]if np.random.random() <= min([Lambda/(state+1), 1]):state = state + 1else:passelse:#选择 n+1 > n，此时接受概率为 min[(n+1)/Lambda, 1]if np.random.random() <= min([(state)/Lambda, 1]):state = state - 1else:passelse:print("undefined state!")breakrandom_list[i] = copy.deepcopy(state)if max_state < state:max_state = copy.deepcopy(state)return random_list, max_statenum = int(1e7)
start = time.time()
results, max_state = poidis(7, num)
end = time.time()
print("average time:", round((end-start)/num, 9))hist_doc = plt.hist(results, bins=[i for i in range(max_state+2)])
print("ave:", np.average(results))
print("std:", np.std(results))
print(hist_doc)plt.show()def block_averages(data, block_size):num_blocks = len(data) // block_sizeblocks = data[:num_blocks*block_size].reshape(num_blocks, block_size)block_avgs = blocks.mean(axis=1)return block_avgsblock_mean = []
block_std  = []for i in range(1, 201):block_size = 5 * iblock_avgs = block_averages(results, block_size)mean_estimate = np.mean(block_avgs)standard_error = np.std(block_avgs, ddof=1) / np.sqrt(len(block_avgs))block_mean.append(mean_estimate)block_std.append(standard_error)plt.scatter(range(1, 201), block_std, s=2)
plt.show()