1

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

1 回答 1

3

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

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

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

于 2013-11-12T13:59:29.807 回答