利用MCMC 获得泊松分布

---

• 写出概率流方程如下

        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()