Esta é a minha primeira tentativa de alguém que vem do campo freqüentador fazer análise de dados bayesiana. Li vários tutoriais e alguns capítulos da análise de dados bayesiana de A. Gelman.

Como o primeiro exemplo de análise de dados mais ou menos independente que escolhi, é o tempo de espera dos trens. Eu me perguntei: qual é a distribuição dos tempos de espera?

O conjunto de dados foi fornecido em um blog e foi analisado de forma ligeiramente diferente e fora do PyMC.

Meu objetivo é estimar os tempos de espera esperados dos trens, considerando essas 19 entradas de dados.

O modelo que eu construí é o seguinte:

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

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

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

$\rho \sim Poisson(\lambda)$

onde é a média dos dados e é o desvio padrão dos dados multiplicado por 1000.$\hat\mu$$\hat\sigma$

Eu modelei o tempo de espera estimado como utilizando a distribuição de Poisson. O parâmetro rate para essa distribuição é modelado usando a distribuição Gamma, pois é uma distribuição conjugada à distribuição Poisson. Os hiperpriores e foram modelados com distribuições Normal e Meia-Normal, respectivamente. O desvio padrão foi feito o mais amplo possível para ser o mais não-compromissivo possível.$\rho$$\mu$$\sigma$$\sigma$

Eu tenho um monte de perguntas

• Esse modelo é razoável para a tarefa (várias maneiras possíveis de modelar?)?
• Cometi algum erro de iniciante?
• O modelo pode ser simplificado (eu tendem a complicar coisas simples)?
• Como posso verificar se o posterior para o parâmetro rate ( ) está realmente ajustando os dados?$\rho$
• Como posso tirar algumas amostras da distribuição de Poisson ajustada para ver as amostras?

As partes posteriores após 5000 etapas da Metropolis são assim:

Também posso postar o código fonte. No estágio de ajuste do modelo, eu faço as etapas para os parâmetros e usando NUTS. Depois, na segunda etapa, faço Metropolis para o parâmetro rate . Finalmente, traço o traço usando as ferramentas incorporadas.$\mu$$\sigma$$\rho$

Ficaria muito grato por quaisquer comentários e comentários que me permitissem entender uma programação mais probabilística. Pode haver exemplos mais clássicos que valem a pena experimentar?

Aqui está o código que escrevi em Python usando PyMC3. O arquivo de dados pode ser encontrado aqui .

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()

Direi primeiro o que faria e depois responderei às perguntas específicas que você tinha.

O que eu faria (pelo menos inicialmente)

Aqui está o que eu recolho do seu post, você tem tempos de espera de treinamento para 19 observações e está interessado em deduzir o tempo de espera esperado.

Vou definir para para ser o tempo de espera para o trem . Não vejo razão para que esses tempos de espera sejam inteiros, portanto, assumirei que sejam quantidades contínuas positivas, ou seja, . Estou assumindo que todos os tempos de espera são realmente observados.$W_i$$i=1,\ldots,19$$i$$W_i\in\mathbb{R}^+$

Existem várias suposições possíveis do modelo que poderiam ser usadas e, com 19 observações, pode ser difícil determinar qual modelo é mais razoável. Alguns exemplos são log-normal, gama, exponencial, Weibull.

Como primeiro modelo, sugiro modelar e, em seguida, assumo Com esta opção, você tem acesso à riqueza da teoria normal que existe, por exemplo, um conjugado anterior. O conjugado anterior é uma distribuição gama-inversa normal, ou seja, que é o distribuição gama inversa. Como alternativa, você pode usar o padrão anterior nesse caso, o posterior também é uma distribuição gama-inversa normal.$Y_i=\log(W_i)$

Como , podemos responder perguntas sobre o tempo de espera esperado, retirando amostras conjuntas de e de sua distribuição posterior, que é inversa normal distribuição de gamma e, em seguida, calculando para cada uma dessas amostras. Isso mostra a partir do posterior para o tempo de espera esperado. μ σ 2 e μ + σ /$E[W_i] = e^{\mu+\sigma^/2}$$\mu$$\sigma^2$$e^{\mu+\sigma^/2}$

Respondendo suas perguntas

• Esse modelo é razoável para a tarefa (várias maneiras possíveis de modelar?)?

Um Poisson não parece apropriado para dados que poderiam ter um valor não inteiro. Você tem apenas um e, portanto, não pode aprender os parâmetros da distribuição gama atribuída a . Outra maneira de dizer isso é que você criou um modelo hierárquico, mas não há estrutura hierárquica nos dados.$\lambda$$\lambda$

• Cometi algum erro de iniciante?

Veja o comentário anterior.

Além disso, será realmente útil se sua matemática e seu código concordarem, por exemplo, onde está nos resultados do MCMC? o que é sd e rate no seu código?$\lambda$

Seu prior não deve depender dos dados.

• O modelo pode ser simplificado (eu tendem a complicar coisas simples)?

Sim e deveria. Veja minha abordagem de modelagem.

• Como posso verificar se o posterior para o parâmetro rate ( ) está realmente ajustando os dados?$\rho$

Não é suposto ser os seus dados? Você quer dizer ? Uma coisa a verificar é garantir que o tempo médio de espera da amostra faça sentido em relação à sua distribuição posterior no tempo médio de espera. A menos que você tenha um histórico bizarro, a média da amostra deve estar próxima do pico da distribuição posterior.$\rho$$\lambda$

• Como posso tirar algumas amostras da distribuição de Poisson ajustada para ver as amostras?

Eu acredito que você quer uma distribuição preditiva posterior. Para cada iteração no seu MCMC, você deve inserir os valores dos parâmetros para essa iteração e coletar uma amostra.