- 写出概率流方程如下
if state == 0:
if np.random.random() <= min([Lambda/2, 1]):
state = 1
else:
pass
elif state == 1:
if choose_prob_state[i] <= 0.5:
#选择 1 -> 0，此时的接受概率为min[2/Lambda, 1]
if np.random.random() <= min([2/Lambda, 1]):
state = 0
else:
pass
else:
#选择 1 -> 2，此时接受概率为 min[Lambda/(n+1), 1]
if np.random.random() <= min([Lambda/(state+1), 1]):
state = 2
else:
pass
elif 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 + 1
else:
pass
else:
#选择 n+1 > n，此时接受概率为 min[(n+1)/Lambda, 1]
if np.random.random() <= min([(state)/Lambda, 1]):
state = state - 1
else:
pass
- blocking 方法
def block_averages(data, block_size):
num_blocks = len(data) // block_size
blocks = data[:num_blocks*block_size].reshape(num_blocks, block_size)
block_avgs = blocks.mean(axis=1)
return block_avgs
block_mean = []
block_std = []
for i in range(1, 201):
block_size = 5 * i
block_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 = init
max_state = init
random_list[0] = state
choose_prob_state = np.random.random(num)
for i in range(1, num):
if state == 0:
if np.random.random() <= min([Lambda/2, 1]):
state = 1
else:
pass
elif state == 1:
if choose_prob_state[i] <= 0.5:
#选择 1 -> 0，此时的接受概率为min[2/Lambda, 1]
if np.random.random() <= min([2/Lambda, 1]):
state = 0
else:
pass
else:
#选择 1 -> 2，此时接受概率为 min[Lambda/(n+1), 1]
if np.random.random() <= min([Lambda/(state+1), 1]):
state = 2
else:
pass
elif 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 + 1
else:
pass
else:
#选择 n+1 > n，此时接受概率为 min[(n+1)/Lambda, 1]
if np.random.random() <= min([(state)/Lambda, 1]):
state = state - 1
else:
pass
else:
print("undefined state!")
break
random_list[i] = copy.deepcopy(state)
if max_state < state:
max_state = copy.deepcopy(state)
return random_list, max_state
num = 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_size
blocks = data[:num_blocks*block_size].reshape(num_blocks, block_size)
block_avgs = blocks.mean(axis=1)
return block_avgs
block_mean = []
block_std = []
for i in range(1, 201):
block_size = 5 * i
block_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()
