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:
gdzie to średnia danych, a to standardowe odchylenie danych pomnożone przez 1000.
Modelowałem oczekiwany czas oczekiwania jako 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 i zostały modelowane odpowiednio z rozkładami Normalny i Pół-Normalny. Sigma odchylenia standardowego 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?
- 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.σ ρ
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()
Odpowiedzi:
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. i = 1 , … , 19 i W i ∈ R +Wi i=1,…,19 i Wi∈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 i n d ∼ N ( μ , σ 2 ) . μ | σ 2 ∼ N ( m , σ 2 C )Yi=log(Wi)
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. μ σ 2 e μ + σ / 2E[Wi]=eμ+σ/2 μ σ2 eμ+σ/2
Odpowiadając na twoje pytania
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.λλ λ
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?λ
Twoje wcześniejsze nie powinny zależeć od danych.
Tak i powinno. Zobacz moje podejście do modelowania.
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.λρ λ
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ę.
źródło