Python中的蒙特卡罗模拟：如何从建模数据中获得完全相同的转移概率？答案

【问题标题】：Monte Carlo simulation in Python: How to get exactly identical transition probabilities from modeled data?Python中的蒙特卡罗模拟：如何从建模数据中获得完全相同的转移概率？
【发布时间】：2013-11-24 03:29:48
【问题描述】：

我已经运行了这个模拟（如下所示）并获得了干到干和湿到湿条件下的模拟转换概率。干到干的模拟结果几乎等于估计的干到干 (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

【问题讨论】：

请不要从您的问题中删除所有代码内容。它使人们难以理解答案所指的内容，并且对于正在寻找类似问题的解决方案的任何人来说都是无用的。

标签： python arrays database for-loop random

【解决方案1】：

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

没有什么可以保证您将获得准确的转换概率 - 在任何一次运行中，您都可能获得任何结果。您所做的是为每个单个转换概率生成一个估计量，其均值等于实际值 (0.345) 和一些正方差。您的估计量的方差随着 n = 样本大小而减小，但始终为正。

如果您想要更接近实际转换概率的值（更快的收敛），请应用一些著名的方差减少技术：Stratified Sampling、Importance Sampling 等 - 太多了，不一一赘述。这是一个快速的技术 - 采用 np.random.rand() 生成的均匀随机偏差，并像往常一样估计。然后使用转换后的偏差生成另一个估计器：[(1-x) for x in stored_deviates]。两个估计量的平均值降低了方差（0.5）。

【讨论】：