To moja pierwsza próba, aby ktoś z obozu dla osób często odwiedzających przeprowadzał analizę danych bayesowskich. Przeczytałem kilka samouczków i kilka rozdziałów z Bayesian Data Analysis autorstwa A. Gelmana.

Jako pierwszy mniej lub bardziej niezależny przykład analizy danych, który wybrałem, są czasy oczekiwania na pociąg. Zadałem sobie pytanie: jaki jest rozkład czasów oczekiwania?

Zestaw danych został udostępniony na blogu i był nieco inaczej analizowany i poza PyMC.

Moim celem jest oszacowanie oczekiwanego czasu oczekiwania na pociąg, biorąc pod uwagę te 19 danych.

Model, który zbudowałem, jest następujący:

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

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

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

$\rho \sim Poisson(\lambda)$

gdzie $\hat\mu$ to średnia danych, a $\hat\sigma$ to standardowe odchylenie danych pomnożone przez 1000.

Modelowałem oczekiwany czas oczekiwania jako $\rho$ przy użyciu rozkładu Poissona. Parametr szybkości dla tego rozkładu jest modelowany przy użyciu rozkładu gamma, ponieważ jest to rozkład sprzężony z rozkładem Poissona. Hiper-priors $\mu$ i $\sigma$ zostały modelowane odpowiednio z rozkładami Normalny i Pół-Normalny. odchylenia standardowego $\sigma$ została tak szeroka, jak to tylko możliwe, aby była możliwie jak najmniej powszechna.

Mam mnóstwo pytań

Czy ten model jest odpowiedni do zadania (kilka możliwych sposobów modelowania?)?
Czy popełniłem jakieś błędy dla początkujących?
Czy można uprościć model (mam tendencję do komplikowania prostych rzeczy)?
Jak mogę sprawdzić, czy tylny parametr szybkości ( ) faktycznie pasuje do danych? $\rho$
Jak narysować próbki z dopasowanego rozkładu Poissona, aby zobaczyć próbki?

Widok boczny po 5000 kroków Metropolis wygląda następująco:

Mogę również opublikować kod źródłowy. Na etapie dopasowania modelu wykonuję kroki dla parametrów i przy użyciu NUTS. Następnie w drugim kroku robię Metropolis dla parametru stawki . Na koniec wykreślam ślad za pomocą wbudowanych narzędzi. $\mu$ $\sigma$ $\rho$

Byłbym bardzo wdzięczny za wszelkie uwagi i komentarze, które umożliwiłyby mi uchwycenie bardziej probabilistycznego programowania. Czy mogą istnieć bardziej klasyczne przykłady, z którymi warto eksperymentować?

Oto kod, który napisałem w Pythonie przy użyciu PyMC3. Plik danych można znaleźć tutaj .

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

bayesian pymc

— Vladislavs Dovgalecs
źródło

Ładne pytanie, ale proponuję edytować tytuł: Twoje pytania są raczej niezależne od oprogramowania i wydają się bardziej na temat oceny modelu. Możesz nawet pokroić go na osobne, powiązane pytania.

— Sean Easter

@SeanEaster Thanks! W rzeczywistości jest to związane z oprogramowaniem, chociaż zgadzam się co do tytułu. Jestem gotów dodać kod źródłowy na żądanie, ponieważ opowiada pełniejszą historię, ale może również sprawić, że pytanie będzie bardziej nieporęczne i potencjalnie bardziej mylące. Edytuj tytuł, ponieważ nie przychodzi mi do głowy nic bardziej ogólnego.

— Vladislavs Dovgalecs

Zgadzam się. Myślę, że to naprawdę dwa pytania. Próbowałem odpowiedzieć na pytania modelujące.

— jaradniemi

Najpierw powiem ci, co bym zrobił, a następnie odpowiem na konkretne pytania, które miałeś.

Co bym zrobił (przynajmniej na początku)

Oto, co zbieram z twojego postu, masz treningowe czasy oczekiwania na 19 obserwacji i jesteś zainteresowany wnioskiem na temat oczekiwanego czasu oczekiwania.

Zdefiniuję dla jako czas oczekiwania na pociąg . Nie widzę powodu, aby te czasy oczekiwania były liczbami całkowitymi, więc że są to dodatnie wielkości ciągłe, tj. . Zakładam, że wszystkie czasy oczekiwania są przestrzegane. $W_i$ $i=1,\ldots,19$ $i$ $W_i\in\mathbb{R}^+$

Istnieje kilka możliwych założeń modelu, które można zastosować, a przy 19 obserwacjach może być trudno ustalić, który model jest bardziej uzasadniony. Niektóre przykłady to log-normal, gamma, wykładniczy, Weibull.

Jako pierwszy model sugerowałbym modelowanie a następnie Dzięki temu wyborowi masz dostęp do bogactwa normalnej teorii, która istnieje, np. Przed koniugatem. Koniugat przed jest rozkładem normalnej odwrotności gamma, tj. gdzie jest odwrotny rozkład gamma. Alternatywnie możesz użyć domyślnego wcześniejszego w którym to przypadku tylny jest również rozkładem normalnej odwrotności gamma. $Y_i=\log(W_i)$

Y_{i} \overset{i n d}{\sim} N (μ, σ^{2}) .

$Y_i \stackrel{ind}{\sim} N(\mu,\sigma^2).$

μ | σ^{2} \sim N (m, σ^{2} C) σ^{2} \sim I G (a, b)

$\mu|\sigma^2 \sim N(m,\sigma^2 C) \quad \sigma^2 \sim IG(a,b)$

I G

$IG$

p (μ, σ^{2}) \propto 1 / σ^{2}

$p(\mu,\sigma^2)\propto 1/\sigma^2$

Ponieważ , możemy odpowiedzieć na pytania dotyczące oczekiwanego czasu oczekiwania, pobierając próbki próbek i z ich rozkładu z tyłu, który jest normalną odwrotnością - rozkład gamma, a następnie obliczenie dla każdej z tych próbek. To próbki z tyłu dla oczekiwanego czasu oczekiwania. $E[W_i] = e^{\mu+\sigma^/2}$ $\mu$ $\sigma^2$ $e^{\mu+\sigma^/2}$

Odpowiadając na twoje pytania

Czy ten model jest odpowiedni do zadania (kilka możliwych sposobów modelowania?)?

Poisson nie wydaje się odpowiedni dla danych, które mogą być wycenione w postaci liczb całkowitych. Masz tylko jedną dlatego nie możesz poznać parametrów rozkładu gamma, które przypisałeś do . Innym sposobem na powiedzenie tego jest zbudowanie modelu hierarchicznego, ale w danych nie ma hierarchicznej struktury. $\lambda$ $\lambda$

Czy popełniłem jakieś błędy dla początkujących?

Zobacz poprzedni komentarz.

Pomoże to również, jeśli matematyka i kod są zgodne, np. Gdzie w wynikach MCMC? co to jest sd i stopa w twoim kodzie? $\lambda$

Twoje wcześniejsze nie powinny zależeć od danych.

Czy można uprościć model (mam tendencję do komplikowania prostych rzeczy)?

Tak i powinno. Zobacz moje podejście do modelowania.

Jak mogę sprawdzić, czy tylny parametr szybkości ( ) faktycznie pasuje do danych? $\rho$

Nie jest ma być dane? Masz na myśli ? Jedną z rzeczy, które należy sprawdzić, jest upewnienie się, że średni czas oczekiwania na próbkę ma sens w stosunku do twojego rozkładu tylnego na średni czas oczekiwania. O ile nie masz dziwnego uprzedniego, średnia próbki powinna znajdować się w pobliżu szczytu rozkładu tylnego. $\rho$ $\lambda$

Jak narysować próbki z dopasowanego rozkładu Poissona, aby zobaczyć próbki?

Wierzę, że chcesz tylnej dystrybucji predykcyjnej. Dla każdej iteracji w MCMC podłączasz wartości parametrów dla tej iteracji i pobierasz próbkę.

— jaradniemi
źródło

Wielkie dzięki! Szybko przeczytałem twoją odpowiedź. Potrzebuję trochę czasu, aby go przetrawić, znaleźć odniesienia do niektórych dystrybucji i pojęć i spróbować zaimplementować go w PyMC. Przy okazji, właśnie dodałem kod Python do mojego eksperymentu.

— Vladislavs Dovgalecs

Bayesowskie modelowanie czasów oczekiwania pociągów: definicja modelu

Co bym zrobił (przynajmniej na początku)

Odpowiadając na twoje pytania