Это моя первая попытка того, чтобы кто-то из лагеря для частых людей сделал анализ байесовских данных. Я прочитал несколько учебников и несколько глав из Байесовского анализа данных А. Гельмана.

В качестве первого более или менее независимого примера анализа данных я выбрал время ожидания поезда. Я спросил себя: каково распределение времени ожидания?

Набор данных был представлен в блоге и анализировался несколько иначе и за пределами PyMC.

Моя цель - оценить ожидаемое время ожидания поезда с учетом этих 19 записей.

Модель, которую я построил, следующая:

$\mu \sim N(\hat\mu,\hat\sigma)$

$\sigma \sim |N(0,\hat\sigma)|$

$\lambda \sim \Gamma(\mu,\sigma)$

$\rho \sim Poisson(\lambda)$

где $\hat\mu$ - среднее значение данных, а $\hat\sigma$ - стандартное отклонение данных, умноженное на 1000.

Я смоделировал ожидаемое время ожидания как $\rho$ используя распределение Пуассона. Параметр скорости для этого распределения моделируется с использованием гамма-распределения, поскольку оно является сопряженным распределением к распределению Пуассона. Гиперприоры $\mu$ и $\sigma$ были смоделированы с нормальным и полунормальным распределениями соответственно. Стандартное отклонение $\sigma$ было сделано как можно более широким, чтобы быть как можно более непринужденным.

У меня есть куча вопросов

• Является ли эта модель приемлемой для данной задачи (несколько возможных способов моделирования?)?
• Я сделал какие-нибудь ошибки новичка?
• Можно ли упростить модель (я склонен усложнять простые вещи)?
• Как я могу проверить, соответствует ли апостериорный параметр скорости ( ) данным?$\rho$
• Как я могу извлечь некоторые образцы из подходящего распределения Пуассона, чтобы увидеть образцы?

Постеры после 5000 шагов Метрополиса выглядят так:

Я также могу опубликовать исходный код. На этапе подбора модели я делаю шаги для параметров и используя NUTS. Затем на втором этапе я делаю Metropolis для параметра скорости . Наконец, я строю трассировку, используя встроенные инструменты.$\mu$$\sigma$$\rho$

Я был бы очень благодарен за любые замечания и комментарии, которые позволили бы мне понять более вероятностное программирование. Может быть, есть более классические примеры, с которыми стоит поэкспериментировать?

Вот код, который я написал на Python, используя PyMC3. Файл данных можно найти здесь .

import matplotlib.pyplot as plt
import pandas as pd
import numpy as np
import pymc3

from scipy import optimize

from pylab import figure, axes, title, show

from pymc3.distributions import Normal, HalfNormal, Poisson, Gamma, Exponential
from pymc3 import find_MAP
from pymc3 import Metropolis, NUTS, sample
from pymc3 import summary, traceplot

df = pd.read_csv( 'train_wait.csv' )

diff_mean = np.mean( df["diff"] )
diff_std = 1000*np.std( df["diff"] )

model = pymc3.Model()

with model:
    # unknown model parameters
    mu = Normal('mu',mu=diff_mean,sd=diff_std)
    sd = HalfNormal('sd',sd=diff_std)

    # unknown model parameter of interest
    rate = Gamma( 'rate', mu=mu, sd=sd )

    # observed
    diff = Poisson( 'diff', rate, observed=df["diff"] )

with model:
    step1 = NUTS([mu,sd])
    step2 = Metropolis([rate])
    trace = sample( 5000, step=[step1,step2] )

plt.figure()
traceplot(trace)
plt.savefig("rate.pdf")
plt.show()
plt.close()

Сначала я расскажу вам, что я буду делать, а затем я отвечу на ваши конкретные вопросы.

Что бы я делал (по крайней мере, изначально)

Вот что я понял из вашего поста: у вас есть время ожидания тренировки для 19 наблюдений, и вы заинтересованы в выводе об ожидаемом времени ожидания.

Я определю для как время ожидания поезда . Я не вижу причин, по которым эти времена ожидания должны быть целыми числами, поэтому я буду предполагать, что они являются положительными непрерывными величинами, то есть . Я предполагаю, что все времена ожидания действительно соблюдаются.$W_i$$i=1,\ldots,19$$i$$W_i\in\mathbb{R}^+$

Существует несколько возможных модельных допущений, которые можно использовать, и с 19 наблюдениями может быть трудно определить, какая модель является более разумной. Некоторые примеры: логарифмическая нормальная, гамма, экспоненциальная, Вейбулла.

В качестве первой модели я бы предложил смоделировать а затем предположить, что С этим выбором у вас есть доступ к богатству обычной теории, которая существует, например, до сопряженного. Сопряженный априор является нормальным гамма-распределением, то есть где - это обратное гамма-распределение. В качестве альтернативы, вы можете использовать предшествующий по умолчанию в этом случае апостериор также является нормальным распределением гамма-инверсии.$Y_i=\log(W_i)$

Поскольку , мы можем ответить на вопросы об ожидаемом времени ожидания, взяв совместные выборки и из их апостериорного распределения, которое является нормальным обратным -гамма-распределение, а затем вычисление для каждого из этих образцов. Это образцы сзади для ожидаемого времени ожидания.$E[W_i] = e^{\mu+\sigma^/2}$$\mu$$\sigma^2$$e^{\mu+\sigma^/2}$

Отвечая на ваши вопросы

• Является ли эта модель приемлемой для данной задачи (несколько возможных способов моделирования?)?

Пуассон не подходит для данных, которые могут быть нецелыми. У вас есть только один и поэтому вы не можете узнать параметры гамма-распределения, которые вы назначили . Другой способ сказать, что вы построили иерархическую модель, но в данных нет иерархической структуры.$\lambda$$\lambda$

• Я сделал какие-нибудь ошибки новичка?

Смотрите предыдущий комментарий.

Кроме того, это действительно поможет, если ваша математика и ваш код согласуются, например, где находится в ваших результатах MCMC? Что такое SD и скорость в вашем коде?$\lambda$

Ваш предыдущий не должен зависеть от данных.

• Можно ли упростить модель (я склонен усложнять простые вещи)?

Да так и должно быть. Смотрите мой подход к моделированию.

• Как я могу проверить, соответствует ли апостериорный параметр скорости ( ) данным?$\rho$

Не должны быть ваши данные? Вы имеете в виду ? Одна вещь, которую нужно проверить, - убедиться, что среднее время ожидания выборки имеет смысл относительно вашего апостериорного распределения среднего времени ожидания. Если у вас нет причудливого предшествующего уровня, среднее значение выборки должно быть близко к пику последующего распределения.$\rho$$\lambda$

• Как я могу извлечь некоторые образцы из подходящего распределения Пуассона, чтобы увидеть образцы?

Я полагаю, что вы хотите последующее прогнозное распределение. Для каждой итерации в MCMC вы вставляете значения параметров для этой итерации и берете образец.