Jak uzyskać prognozę dla określonej zmiennej w WinBUGS?

10

Jestem nowym użytkownikiem WinBUGS i mam jedno pytanie do twojej pomocy. Po uruchomieniu następującego kodu uzyskałem parametry beta0through beta4(statystyki, gęstość), ale nie wiem, jak uzyskać prognozę ostatniej wartości h, którą ustawiłem NAdo modelowania w kodzie.

Czy ktoś może mi podpowiedzieć? Wszelkie porady będą mile widziane.


model {
for(i in 1: N) {
CF01[i] ~ dnorm(0, 20)
CF02[i]  ~ dnorm(0, 1)
h[i] ~ dpois (lambda [i])
log(lambda [i]) <- beta0 + beta1*CF03[i] + beta2*CF02[i] + beta3*CF01[i] + beta4*IND[i]
}
beta0 ~ dnorm(0.0, 1.0E-6)
beta1 ~ dnorm(0.0, 1.0E-6)
beta2 ~ dnorm(0.0, 1.0E-6)
beta3 ~ dnorm(0.0, 1.0E-6)
beta4  <- log(p)
p ~ dunif(lower, upper)
}

INITS
list(beta0 = 0, beta1 = 0, beta2 = 0, beta3 = 0, p = 0.9)

DATA(LIST)
list(N = 154, lower = 0.80, upper = 0.95,

h = c(1, 4, 1, 2, 1, 2, 1, 1, 1, 3, 3, 0, 0, 0, 2, 0, 1, 0, 4, 2,
3, 0, 2, 1, 1, 2, 2, 2, 3, 4, 2, 3, 1, 0, 1, 3, 3, 3, 1, 0, 1,
0, 5, 2, 1, 2, 1, 3, 3, 1, 1, 0, 2, 2, 0, 3, 0, 0, 3, 2, 2, 2,
1, 0, 3, 3, 1, 1, 1, 2, 1, 0, 1, 2, 1, 2, 0, 2, 1, 0, 0, 2, 5,
0, 2, 1, 0, 2, 1, 2, 2, 2, 0, 3, 2, 1, 3, 3, 3, 3, 0, 1, 3, 3,
3, 1, 0, 0, 1, 2, 1, 0, 1, 4, 1, 1, 1, 1, 2, 1, 3, 0, 0, 1, 1,
1, 1, 0, 2, 1, 0, 0, 1, 1, 5, 1, 1, 1, 3, 0, 1, 1, 1, 0, 2, 1,
0, 3, 3, 0, 0, 1, 2, 6, NA),

CF03 = c(-1.575, 0.170, -1.040, -0.010, -0.750,
0.665, -0.250, 0.145, -0.345, -1.915, -1.515,
0.215, -1.040, -0.035, 0.805, -0.860, -1.775,
1.725, -1.345, 1.055, -1.935, -0.160, -0.075,
-1.305, 1.175, 0.130, -1.025, -0.630, 0.065,
-0.665, 0.415, -0.660, -1.145, 0.165, 0.955,
-0.920, 0.250, -0.365, 0.750, 0.045, -2.760,
-0.520, -0.095, 0.700, 0.155, -0.580, -0.970,
-0.685, -0.640, -0.900, -0.250, -1.355, -1.330,
0.440, -1.505, -1.715, -0.330, 1.375, -1.135,
-1.285, 0.605, 0.360, 0.705, 1.380, -2.385, -1.875,
-0.390, 0.770, 1.605, -0.430, -1.120, 1.575, 0.440,
-1.320, -0.540, -1.490, -1.815, -2.395, 0.305,
0.735, -0.790, -1.070, -1.085, -0.540, -0.935,
-0.790, 1.400, 0.310, -1.150, -0.725, -0.150,
-0.640, 2.040, -1.180, -0.235, -0.070, -0.500,
-0.750, -1.450, -0.235, -1.635, -0.460, -1.855,
-0.925, 0.075, 2.900, -0.820, -0.170, -0.355,
-0.170, 0.595, 0.655, 0.070, 0.330, 0.395, 1.165,
0.750, -0.275, -0.700, 0.880, -0.970, 1.155, 0.600,
-0.075, -1.120, 1.480, -1.255, 0.255, 0.725,
-1.230, -0.760, -0.380, -0.015, -1.005, -1.605,
0.435, -0.695, -1.995, 0.315, -0.385, -0.175,
-0.470, -1.215, 0.780, -1.860, -0.035, -2.700,
-1.055, 1.210, 0.600, -0.710, 0.425, 0.155, -0.525,
-0.565),

CF02 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, 0.38, 0.06, -0.94,
-0.02, -0.28, -0.78, -0.95, 2.33, 1.43, 1.24, 1.26,
-0.75, -1.5, -2.09, 1.01, -0.05, 2.48, 2.48, 0.46,
0.46, -0.2, -1.11, 0.52, -0.37, 0.58, 0.86, 0.59,
-0.12, -1.33, 1.4, -1.84, -1.4, -0.76, -0.23,
-1.78, -1.43, 1.2, 0.32, 1.87, 0.43, -1.71, -0.54,
-1.25, -1.01, -1.98, 0.52, -1.07, -0.44, -0.24,
-1.31, -2.14, -0.43, 2.47, -0.09, -1.32, -0.3,
-0.99, 1.1, 0.41, 1.01, -0.19, 0.45, -0.07, -1.41,
0.87, 0.68, 1.61, 0.36, -1.06, -0.44, -0.16, 0.72,
-0.69, -0.94, 0.11, 1.25, 0.33, -0.05, 0.87, -0.37,
-0.2, -2.22, 0.26, -0.53, -1.59, 0.04, 0.16, -2.66,
-0.21, -0.92, 0.25, -1.36, -1.62, 0.61, -0.2, 0,
1.14, 0.27, -0.64, 2.29, -0.56, -0.59, 0.44, -0.05,
0.56, 0.71, 0.32, -0.38, 0.01, -1.62, 1.74, 0.27, 0.97,
1.22, -0.21, -0.05, 1.15, 1.49, -0.15, 0.05, -0.87,
-0.3, -0.08, 0.5, 0.84, -1.67, 0.69, 0.47, 0.44,
-1.35, -0.24, -1.5, -1.32, -0.08, 0.76, -0.57,
-0.84, -1.11, 1.94, -0.68),

CF01 = c(NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, NA, NA, NA, NA, NA, NA, NA, NA, NA,
NA, -0.117, -0.211, -0.333, -0.229, -0.272,
-0.243, -0.148, 0.191, -0.263, -0.239, -0.168,
-0.381, -0.512, -0.338, -0.296, 0.067, 0.104,
-0.254, -0.167, -0.526, -0.096, -0.43, 0.013,
-0.438, -0.297, -0.131, -0.098, -0.046, -0.063,
-0.194, -0.155, -0.645, -0.603, -0.374, -0.214,
-0.165, -0.509, -0.171, -0.442, -0.468, -0.289,
-0.427, -0.519, -0.454, 0.046, -0.275, -0.401,
-0.542, -0.488, -0.52, -0.018, -0.551, -0.444,
-0.254, -0.286, 0.048, -0.03, -0.015, -0.219,
-0.029, 0.059, 0.007, 0.157, 0.141, -0.035, 0.136,
0.526, 0.113, 0.22, -0.022, -0.173, 0.021, -0.027,
0.261, 0.082, -0.266, -0.284, -0.097, 0.097, -0.06,
0.397, 0.315, 0.302, -0.026, 0.268, -0.111, 0.084,
0.14, -0.073, 0.287, 0.061, 0.035, -0.022, -0.091,
-0.22, -0.021, -0.17, -0.184, 0.121, -0.192,
-0.24, -0.283, -0.003, -0.45, -0.138, -0.143,
0.017, -0.245, 0.003, 0.108, 0.015, -0.219, 0.09,
-0.22, -0.004, -0.178, 0.396, 0.204, 0.342, 0.079,
-0.034, -0.122, -0.24, -0.125, 0.382, 0.072, 0.294,
0.577, 0.4, 0.213, 0.359, 0.074, 0.388, 0.253, 0.167),

IND = c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1,
1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
0, 0, 0, 0, 0, 0, 0, 0, 0, 0))
Bo Yu
źródło
Czy nie pytasz tylko o wartość lambda [N]?
whuber
@ whuber tak, myślę, że to prawda, ale bardziej fundamentalnie, musisz mieć rzeczy, które chcesz przewidzieć (tj. mieć rozkład późniejszy), odróżniać się od rzeczy, które już zaobserwowałeś. Możesz dokonać prognozy wprost w winbugs lub w postprocessingu, korzystając z próbek beta.
atiretoo
@atiretoo O ile mi wiadomo, lambdy są dokładnie tym, co chce się przewidzieć: jest to uogólniony model liniowy dla rozkładu Poissona z łączem log, a lambda są przewidywanymi parametrami Poissona. Nie zostały zaobserwowane. Uważam, że wszystko, co musisz tutaj zrobić, to ustawić monitor na lambda [N].
whuber
@ Whuber, wolę powiedzieć monitor h[N]zamiast lambda[N]... a dostaniesz rozkład tylny przewidywanej wartości.
Ciekawy
@ tomek, ale h[N]nie jest przewidywaną wartością: będzie to zbiór losowań z zestawu przewidywanych rozkładów Poissona. Jako taki łączy zmienność parametrów Poissona i odchylenie od samych rozkładów Poissona. Istotny jest rozkład tylny lambda[N].
whuber

Odpowiedzi:

6

Wystarczy dodać zmienną hdo listy parametrów, które mają być monitorowane. Jeśli używasz pakietu takiego jak R2WinBUGS, dodaj zmienną hdo listy przekazywanej do parameters.to.saveargumentu bugsfunkcji. Następnie spójrz na swoją ostatnią wartość w h(tą z NA) - dostaniesz tam rozkład boczny.

Jest to zwykły sposób dokonywania prognoz na podstawie wnioskowania bayesowskiego ( patrz także to pytanie ). To jest miłe i proste! Nigdy więcej rozdzielania oceny parametrów i prognoz. Wszystko odbywa się na raz. Późniejszy rozkład parametrów jest podawany przez rzeczywiste dane i propagowany do wartości NA (jako „prognozy”).

Ciekawy
źródło
Tomas, dzięki za pomoc. Próbuję monitorować zmienną hw narzędziu Sample Monitor Tool, ale to nie działa. Czy możesz mi jeszcze raz pomóc? Oto procedura, którą wykonałem w WinBUGS (nie wiem, jak korzystać z R2WinBUGS): 1) wybierz Próbkę w narzędziu Sample Monitor Tool 2) wpisz h w białym polu oznaczonym węzłem 3) kliknij przycisk oznaczony zestaw 4) h jest nie na liście parametrów, które chcę monitorować, podczas gdy inne parametry (beta0, beta1, beta2, beta3, p) są wyświetlane na liście. Czy wiesz, jak mogę dodać „h” do listy parametrów, które chcę monitorować? Dzięki raz jeszcze!
Bo Yu
@ Boou, nie wiem jak to zrobić bezpośrednio w WinBUGS, ponieważ uruchamiam WinBUGS od R, używając pakietu R2WinBUGS. Jest to o wiele bardziej praktyczne, ponieważ możesz po prostu zapisać skrypt R i uruchomić go jako partię, wraz z tworzeniem własnych wykresów itp. Spójrz tutaj na przykład skrypty.
Ciekawy
To powiedziawszy, z pewnością będzie to również możliwe w samym WinBUGS, ale nie wiem jak (i ​​chyba większość ludzi nazywa to od R).
Ciekawy
1
Przede wszystkim dziękuję, whuber, atiretoo i Tomas! Jak już wspomniano wcześniej, tak, jest to uogólniony model liniowy, zmienna h jest dopasowywana przez rozkład Poissona ze zmienną szybkością (lambda) uwarunkowaną różnymi predyktorami (CF01, CF02, CF03 i IND). Ostatnia wartość h jest tym, co muszę wiedzieć i nie jest przestrzegana (oznaczona jako NA), podczas gdy wszystkie inne wartości h są obserwowane. Myślę, że whuber ma rację, muszę ustawić lambda jako parametr w narzędziu Sample Monitor Tool i sprawdzić statystyki ostatniej wartości lambda, a następnie uzyskać moje przewidywania dotyczące ostatniej h. Dziękuje wszystkim.
Bo Yu
1
@ Tomas, dziękuję bardzo. Tak masz rację! WinBUGS zapewnia prognozę h [N], w tym statystyki i gęstość prawdopodobieństwa. Mam to teraz. Z pozdrowieniami,
Bo Yu,