Python中的montecarlo模拟：如何从建模数据中获得完全相同的转移概率？

Question

我已经运行了这个模拟（如下所示），得到了干燥到干燥和潮湿到潮湿条件下的模拟转换概率。干转干的模拟结果几乎与估计的干到干（d2d_tran）相等。但是，模拟的湿到湿值远远低于估计值。程序好像有问题。我试过其他几种方法，但都没有达到预期的效果。你能运行程序并建议我如何提高湿到湿概率的结果吗？提前谢谢。在

我的代码：

import numpy as np
import random, datetime

d2d = np.zeros(12)
d2w = np.zeros(12)
w2w = np.zeros(12)
w2d = np.zeros(12)
pd2d = np.zeros(12)
pw2w = np.zeros(12)

dry = [0.333] ##unconditional probability of dry for January
d2d_tran = [0.564,0.503,0.582,0.621,0.634,0.679,0.738,0.667,0.604,0.564,0.577,0.621]
w2w_tran = [0.784,0.807,0.8,0.732,0.727,0.728,0.64,0.64,0.665,0.717,0.741,0.769]
mu = [3.71,4.46,4.11,2.94,3.01,2.87,2.31,2.44,2.56,3.45,4.32,4.12]
sigma = [6.72,7.92,7.49,6.57,6.09,5.53,4.38,4.69,4.31,5.71,7.64,7.54]

days = np.array([31,28,31,30,31,30,31,31,30,31,30,31])
rain = np.array([])

for y in xrange(0,10000):
    for m in xrange(0,12):
    #Include leap years in the calculation and creat random variables for each month
        if ((y%4 == 0 and y%100 != 0) or y%400 == 0) and m==1: 
            random_num = np.random.rand(29)
        else:
            random_num = np.random.rand(days[m])

    #lets generate a rainfall amount for first day of the random series
        if random_num[0] <= dry[0]:
            random_num[0] = 0
        else:
            random_num[0] = abs(random.gauss(mu[0],sigma[0]))

# generate the whole series in sequence of month and year
        for i in xrange(0,days[m]):
            if random_num[i-1] == 0: #if yesterday was dry
                if random_num[i] <= d2d_tran[m]: #check today against the dry2dry transition probabilities
                    random_num[i] = 0
                    d2d[m] += 1.0
                else:
                    random_num[i] = abs(random.gauss(mu[m],sigma[m]))
                    d2w[m] += 1.0

            else:
                if random_num[i] <= w2w_tran[m]:
                    random_num[i] = abs(random.gauss(mu[m],sigma[m]))
                    w2w[m] += 1.0                         
                else:
                    random_num[i] = 0
                    w2d[m] += 1.0



        pd2d[m] = d2d[m]/(d2d[m] + d2w[m])
        pw2w[m] = w2w[m]/(w2d[m] + w2w[m])


print 'Simulated transition probability of dry2dry:\n', np.around(pd2d, decimals=3)
print 'Simulated transition probability of wet2wet:\n', np.around(pw2w, decimals=3)

### pd2d and pw2w of generated data should be identical to d2d_tran and w2w_tran respectively

Answer 1

这个模拟看起来是正确的，在运行了8000年之后，大多数情况下，我得到的转换概率都在0.001以内，并且随着天数的增加，会有收敛性。在

没有什么能保证你会得到确切的转换概率-在任何一次运行中你都可能得到任何东西。您所做的是为每个单一的转移概率生成一个估计量，它的平均值等于实际值（0.345），以及一些正的方差。估计量的方差随着n=样本量的增加而减小，但它始终是正的。在

如果您希望值更接近实际的转换概率（更快的收敛），请应用一些众所周知的方差缩减技术：Stratified Sampling，Importance Sampling等等，太多了，不必赘述。这里有一个快速的技术-取由np.随机.rand（），照常估算。然后使用转换后的偏差生成另一个估计器：[(1-x) for x in stored_deviates]。两个估计量的平均值使方差减少了0.5。在