Как получить прогноз для конкретной переменной в WinBUGS?

10

Я новый пользователь WinBUGS и у меня есть один вопрос для вашей помощи. После запуска следующего кода я получил параметры beta0сквозного beta4(статистика, плотность), но я не знаю, как получить прогноз последнего значения h, которое я установил NAдля моделирования в коде.

Кто-нибудь может дать мне подсказку? Любой совет будет принята с благодарностью.

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))
bayesian bugs prediction winbugs Бо Ю
источник
Разве вы не спрашиваете о значении лямбды [N]?
whuber
@whuber да, я думаю, что это правильно, но, что более важно, вам нужно, чтобы вещи, которые вы хотите предсказать (т.е. иметь апостериорное распределение), отличались от вещей, которые вы уже наблюдали. Вы можете сделать прогноз явно в winbugs или в постобработке, используя образцы бета-версий.
atiretoo - восстановить монику
@atiretoo Насколько я могу судить, лямбды - это именно то, что каждый хочет предсказать: это обобщенная линейная модель для распределения Пуассона с логарифмической связью, а лямбды - это предсказанные параметры Пуассона. Они не были замечены. Я считаю, что все, что нужно сделать здесь, это установить монитор на лямбду [N].
whuber
@whuber, я бы лучше сказал монитор h[N]вместо lambda[N]... и вы получили апостериорное распределение прогнозируемого значения.
любопытно
@tomek, но h[N]не является прогнозируемым значением: это будет набор отрисовок из набора прогнозируемых распределений Пуассона. Как таковой он сочетает в себе изменение параметров Пуассона и изменение самих этих распределений Пуассона. Что важно, так это заднее распределение lambda[N].
whuber

Ответы:

6

Просто добавьте переменную hв список контролируемых параметров. Если вы используете такой пакет, как R2WinBUGS, добавьте переменную hв список, передаваемый parameters.to.saveаргументу bugsфункции. Затем посмотрите на ваше последнее значение в h(тот, что с NA) - вы получите апостериорное распределение там.

Это обычный способ делать прогнозы в байесовском выводе ( см. Также этот вопрос ). Это красиво и просто! Нет больше разделения оценки параметров и прогнозирования. Все сделано сразу. Последующее распределение параметров задается фактическими данными и распространяется на значения NA (как «прогнозы»).

любознательный
источник
Томас, спасибо за вашу помощь. Я пытаюсь контролировать переменную h в Sample Monitor Tool, но она не работает. Не могли бы вы помочь мне снова? Ниже описана процедура, которую я сделал в WinBUGS (я не знаю, как использовать R2WinBUGS): 1) выберите Sample in Sample Monitor Tool 2) введите h в белом поле, помеченном узлом 3) нажмите на кнопку с пометкой set 4) h is нет в списке параметров, которые я хочу отслеживать, в то время как другие параметры (бета0, бета1, бета2, бета3, р) отображаются в списке. Вы знаете, как я могу добавить «h» в список параметров, которые я хочу контролировать? Еще раз спасибо
Бо Ю
@ BoYu, я не знаю, как сделать это напрямую в WinBUGS, так как я запускаю WinBUGS из R, используя пакет R2WinBUGS. Это гораздо более практично, потому что вы можете просто сохранить сценарий R и запустить его как пакет, наряду с созданием собственных графиков и т. Д. Посмотрите здесь примеры сценариев.
любопытно
Тем не менее, это наверняка будет возможно и в самой WinBUGS, но я не знаю как (и я думаю, что большинство людей называют это из R).
любопытно
1
Прежде всего, спасибо тебе, whuber, atiretoo и Tomas! Как уже упоминалось, да, это обобщенная линейная модель, переменная h подгоняется распределением Пуассона с переменной скоростью (лямбда), обусловленной различными предикторами (CF01, CF02, CF03 и IND). Последнее значение h - это то, что мне нужно знать, и оно не наблюдается (помечено как NA), тогда как все остальные значения h соблюдаются. Я думаю, что правильно, мне нужно установить лямбду в качестве параметра в Sample Monitor Tool и проверить статистику последнего значения лямбды, а затем получить то, что я предсказываю за последние h. Спасибо всем.
Бо Ю
1
@ Томас, спасибо тебе большое. Да ты прав! WinBUGS обеспечивает прогноз h [N], включая статистику и плотность вероятности. Я понял. С наилучшими пожеланиями,
Бо Юй