Последовательный Монте-Карло

Мне дали эту модель, и чтобы получить вероятность, я должен смоделировать данные.

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 следует нормальному распределению.

Я попытался инвертировать функцию, но не могу этого сделать, потому что понятия не имею, будет ли мой крестик положительным или отрицательным. Я понимал, что должен использовать последовательный монте-карло, но не могу понять, как найти функции алгоритма. Что такое 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 montecarlo distribution normal-distribution
person Frank    schedule 15.11.2017    source источник

Ответы (1)


arrow_upward
1
arrow_downward

С такой моделью временных рядов, как ваша, по сути, единственный способ вычислить распределение вероятностей 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()

что дает: введите описание изображения здесь

person rwp    schedule 17.02.2018