Извлечение коэффициентов из объекта mcmc - PullRequest
Новая
       15

Извлечение коэффициентов из объекта mcmc

0 голосов
/ 05 ноября 2018

Извлечение материала из объектов всегда было одним из самых запутанных аспектов R для меня. Я установил байесовскую модель линейной регрессии, используя rjags, и у меня есть следующий объект mcmc: 

summary(m_csim)
Iterations = 1:150000
Thinning interval = 1 
Number of chains = 1 
Sample size per chain = 150000 

1. Empirical mean and standard deviation for each variable,
   plus standard error of the mean:

            Mean        SD  Naive SE Time-series SE
BR2     0.995805 0.0007474 1.930e-06      3.527e-06
BR2adj  0.995680 0.0007697 1.987e-06      3.633e-06
b[1]   -5.890842 0.1654755 4.273e-04      1.289e-02
b[2]    1.941420 0.0390239 1.008e-04      1.991e-03
b[3]    1.056599 0.0555885 1.435e-04      5.599e-03
sig2    0.004678 0.0008333 2.152e-06      3.933e-06

2. Quantiles for each variable:

            2.5%       25%       50%       75%    97.5%
BR2     0.994108  0.995365  0.995888  0.996339  0.99702
BR2adj  0.993932  0.995227  0.995765  0.996229  0.99693
b[1]   -6.210425 -6.000299 -5.894810 -5.784082 -5.55138
b[2]    1.867453  1.914485  1.940372  1.967466  2.02041
b[3]    0.942107  1.020846  1.057720  1.094442  1.16385
sig2    0.003321  0.004082  0.004585  0.005168  0.00657

Чтобы извлечь коэффициенты, я сделал b = colMeans(mod_csim)[3:5]. Я хочу рассчитать достоверные интервалы, поэтому мне нужно также извлечь квантили 0,025 и 0,975. Как мне сделать это программно?

Ответы [ 2 ]

0 голосов
/ 07 ноября 2018

Вероятно, вы можете извлечь квантили напрямую. Как уже отмечали другие, вы можете вызвать str(m_csim), и вы можете ограничить вывод вызова str с помощью str(m_csim, max.level=1) и продолжать добавлять его к аргументу max.level=, пока не увидите нечто, похожее на квантили.

Мне нравится преобразовывать вывод MCMC в data.frame, чтобы с ним было проще работать. Я использую jagsUI, а не rjags, но я часто делаю что-то вроде 

mcmc_df <- as.data.frame(as.matrix(MY_MCMC_OBJECT$samples))

Примечание: это может быть немного по-другому с rjags, но я уверен, что вы можете найти это, немного покопавшись.

Преимущество: я могу получить доступ к одному вектору для каждого mcmc_df$PARAMETER и создать матрицу квантилей с 

mcmc_quants <- apply(mcmc_df, 2, quantile, p=c(0.025, 0.25, 0.5, 0.75, 0.975))

или любые другие квантили, которые вы хотите.

0 голосов
/ 05 ноября 2018

Я надеюсь, что я не переступаю границы своих знаний, но я хочу ответить с точки зрения «в целом», а не специально для рьягов. m_csim является объектом, и многие методы могут быть использованы на нем. Вы использовали метод резюме, чтобы увидеть что-то. Как прокомментировали люди, вероятно, существует метод коэффи. Однако, как прокомментировал кто-то другой (пока я отвечал!), Использование str (), чтобы увидеть, что содержится в объекте, является лучшим способом узнать, какая информация находится в объекте и как к ней обращаться. Я был бы очень удивлен, если бы использование str () не показывало, как найти не только коэффициенты, но и достаточно информации о доверительных интервалах, чтобы позволить вам найти желаемый CI.

...