Установите байесовскую линейную регрессию и предскажите ненаблюдаемые значения.

Question

Я хотел бы использовать Jags plus R для настройки линейной модели с наблюдаемыми величинами и сделать вывод о ненаблюдаемых. Я нашел много примеров в Интернете о том, как настроить модель, но ничего о том, как экстраполировать ее коэффициенты после установки модели в среде Jags. Поэтому я буду признателен за любую помощь в этом.

Мои данные выглядит следующим образом:

ngroups <- 2 
group <- 1:ngroups 
nobs <- 100 
dta <- data.frame(group=rep(group,each=nobs),y=rnorm(nobs*ngroups),x=runif(nobs*ngroups)) 
head(dta) 

Answer 1

JAGS имеет мощные способы, чтобы сделать вывод о недостающих данных, и как только вы получите повесить ее, это просто! Я настоятельно рекомендую вам ознакомиться с 4444975335667294988 Marc Kéry, который представляет собой замечательное введение в программирование на языке BUGS (JAGS достаточно близко к BUGS, что почти все передает).

Простейший способ сделать это включает, как вы говорите, настройку модели. Ниже я предоставляю полный рабочий пример того, как это работает. Но вы, кажется, просите способ получить интервал предсказания без повторного запуска модели (ваша модель очень большая и вычислительно дорогая?). Это также можно сделать.
Как предсказать - трудный путь (без повторного запуска модели) Для каждой итерации MCMC имитируйте отклик для желаемого значения x, основанного на задних рисунках этой итерации для значений ковариации. Итак, представьте, что вы хотите предсказать значение для X = 10. Тогда, если итерации 1 (после приработки) имеет наклон = 2, перехватывать = 1, а стандартное отклонение = 0,5, нарисовать Y-значение из

Y=rnorm(1, 1+2*10, 0.5) 

и повторить для итерации 2, 3, 4, 5. .. Это ваши задние розыгрыши для ответа при X = 10. Примечание:, если вы не отслеживали стандартное отклонение в модели JAGS, вам не повезло и вам нужно снова подгонять модель.

Как предсказать - легкий путь - с работавшего например Основная идея заключается во вставке (в данных) х-значений, чьи ответы вы хотите, чтобы предсказать, с соответствующими у значений NA. Например, если вам нужен интервал прогнозирования для X = 10, вам просто нужно включить точку (10, NA) в свои данные и установить монитор трассировки для значения y.

Я использую JAGS из R с пакетом rjags. Ниже приведен полный рабочий пример, который начинается с моделирования данных, затем добавляет некоторые дополнительные значения х в данные, задает и запускает линейную модель в JAGS через rjags и суммирует результаты. Y [101: 105] содержит рисунки с интервалов заднего предсказания для X [101: 105]. Обратите внимание: Y [1: 100] просто содержит значения y для X [1: 100]. Это наблюдаемые данные, которые мы подавали в модель, и они никогда не меняются по мере обновления модели.

library(rjags) 
# Simulate data (100 observations) 
my.data <- as.data.frame(matrix(data=NA, nrow=100, ncol=2)) 
names(my.data) <- c("X", "Y") 
# the linear model will predict Y based on the covariate X 

my.data$X <- runif(100) # values for the covariate 
int <- 2  # specify the true intercept 
slope <- 1 # specify the true slope 
sigma <- .5 # specify the true residual standard deviation 
my.data$Y <- rnorm(100, slope*my.data$X+int, sigma) # Simulate the data 

#### Extra data for prediction of unknown Y-values from known X-values 
y.predict <- as.data.frame(matrix(data=NA, nrow=5, ncol=2)) 
names(y.predict) <- c("X", "Y") 
y.predict$X <- c(-1, 0, 1.3, 2, 7) 

mydata <- rbind(my.data, y.predict) 


set.seed(333) 
setwd(INSERT YOUR WORKING DIRECTORY HERE) 
sink("mymodel.txt") 
cat("model{ 

    # Priors 

    int ~ dnorm(0, .001) 
    slope ~ dnorm(0, .001) 
    tau <- 1/(sigma * sigma) 
    sigma ~ dunif(0,10) 

    # Model structure 

    for(i in 1:R){ 
    Y[i] ~ dnorm(m[i],tau) 
    m[i] <- int + slope * X[i] 
    } 
    }", fill=TRUE) 
sink() 
jags.data <- list(R=dim(mydata)[1], X=mydata$X, Y=mydata$Y) 

inits <- function(){list(int=rnorm(1, 0, 5), slope=rnorm(1,0,5), 
         sigma=runif(1,0,10))} 

params <- c("Y", "int", "slope", "sigma") 

nc <- 3 
n.adapt <-1000 
n.burn <- 1000 
n.iter <- 10000 
thin <- 10 
my.model <- jags.model('mymodel.txt', data = jags.data, inits=inits, n.chains=nc, n.adapt=n.adapt) 
update(my.model, n.burn) 
my.model_samples <- coda.samples(my.model,params,n.iter=n.iter, thin=thin) 
summary(my.model_samples)