顺序蒙特卡罗答案

【问题标题】：Sequential Monte Carlo顺序蒙特卡罗
【发布时间】：2018-04-28 12:05:44
【问题描述】：

我得到了这个模型，并获得了我应该模拟数据的概率。

x_1 ∼N(0, 102)

x_t =0.5 ∗ (x_t−1) + 25 · (x_t−1)/(1 + (x_t-1)^2) + 8 · cos(1.2 ∗ (t − 1)) + εt
, t = 2, 3, ..

y_t =(x_t)^2/25 + ηt, t = 1, 2, 3, ...

其中 εT 和 ηt 服从正态分布。

我试图反转函数，但我不能这样做，因为我不知道我的 X 是正数还是负数。我知道我应该使用顺序蒙特卡罗，但我不知道如何找到算法的功能。 f 和 g 是什么，我们如何确定 x(t-1) 是否因为 x 的平方而同样可能是正数或负数？

算法：

1 Sample X1 ∼ g1(·). Let w1 = u1 = f1(x1)/g1(x1). Set t = 2

2 Sample Xt|xt−1 ∼ gt(xt|xt−1).

3 Append xt to x1:t−1, obtaining xt

4 Let ut = ft(xt|xt−1)/gt(xt|xt−1)

5 Let wt = wt−1ut , the importance weight for x1:t

6 Increment t and return to step 2

【问题讨论】：

标签： r statistics distribution normal-distribution montecarlo

【解决方案1】：

对于像你这样的时间序列模型，基本上计算 x 或 y 概率分布的唯一方法是运行模型的多次模拟，随机绘制 x_0、eps_t、eta_t 值，然后通过以下方式构建直方图聚合所有运行的样本。在非常特殊的情况下（例如阻尼布朗运动），可能可以通过代数方式计算得到的概率分布，但我认为您的模型没有任何机会。

在 Python 中（恐怕我对 R 语言不够流利），您可以通过以下方式模拟时间序列：

import math, random    

def simSequence(steps, eps=0.1, eta=0.1):
    x = random.normalvariate(0, 102)

    ySamples = []

    for t in range(steps):
        y = (x ** 2) / 25 + random.normalvariate(0, eta)
        ySamples.append(y)

        x = (0.5 * x + 25 * x / (1 + x ** 2)
                + 8 * math.cos(1.2 * t) + random.normalvariate(0, eps))

    return ySamples

（这会将您的 t=1..n 替换为 t=0..(n-1)。）

然后您可以生成 y 时间序列的几个示例的图：

import matplotlib.pyplot as plt

nSteps = 100
for run in range(5):
    history = simSequence(nSteps)
    plt.plot(range(nSteps), history)
plt.show()

得到类似的东西：

如果您想计算 y 在不同时间的概率分布，您可以生成一个矩阵，其列表示 y_t 在公共时间值的实现，并在选定的 t 值处计算直方图：

import numpy

runs = numpy.array([ simSequence(nSteps) for run in range(10000) ])
plt.hist(runs[:,5], bins=25, label='t=5', alpha=0.5, normed=True)
plt.hist(runs[:,10], bins=25, label='t=10', alpha=0.5, normed=True)
plt.legend(loc='best')
plt.show()

给出：

【讨论】：