Трансформация переменных в рдане (байесовский анализ)

Question

Я новичок в байесовском анализе и стараюсь использовать rstan для оценки распределения плотности задних частот. Упражнение пытается воссоздать пример, предоставленный нам моим университетом, используя stan, но я немного смущен относительно того, как правильно преобразовывать переменные. Мой текущий код работает без ошибок, но результат не соответствует тому, что мы предоставили университету (хотя и близко), ниже на графике ниже для ясности с оценками stan в черном. Я получил код для работы, посоветовавшись с руководством и собирая случайные биты, но, в частности, я не уверен, почему нужен target, и если преобразование гаммы действительно верно. Любое руководство будет оценено!

Модель

Стан Код

data { 
    int<lower=0> I; 
    int<lower=0> n[I]; 
    int<lower=0> x[I]; 
    real<lower=0> a; 
    real<lower=0> b; 
    real m; 
    real<lower=0> p; 
} 

parameters { 
    real<lower=0> lambda; 
    real mu; 
    real<lower=0, upper=1> theta[I]; 
} 

transformed parameters { 
    real gam[I]; 
    for(j in 1:I) 
    gam[j] = log(theta[j]/(1-theta[j])) ; 
} 


model { 
    target += gamma_lpdf(lambda | a, b); 
    target += normal_lpdf(mu | m , 1/sqrt(p)); 
    target += normal_lpdf(gam | mu, 1/sqrt(lambda)); 
    target += binomial_lpmf(x | n , theta); 
} 

R код

library(rstan) 
fit <- stan(
    file = "hospital.stan" , 
    data = dat , 
    iter = 20000, 
    warmup = 2000, 
    chains = 1 
) 

Дат

structure(
    list(
     I = 12L, 
     n = c(47, 211, 810, 148, 196, 360, 119, 207, 97, 256, 148, 215), 
     x = c(0, 8, 46, 9, 13, 24, 8, 14, 8, 29, 18, 31), 
     a = 2, 
     b = 2, 
     m = 0, 
     p = 0.01), 
    .Names = c("I", "n", "x", "a", "b", "m", "p") 
) 

--- UPDATE раствором ---

Этот вопрос, как отметил Бен Goodrich, что я получения гамма от тета, где, как это должно было быть другим так как гамма - это моя случайная величина. Правильный строковый код приведен ниже.

data { 
    int<lower=0> I; 
    int<lower=0> n[I]; 
    int<lower=0> x[I]; 
    real<lower=0> a; 
    real<lower=0> b; 
    real m; 
    real<lower=0> p; 
} 

parameters { 
    real<lower=0> lambda; 
    real mu; 
    real gam[I]; 
} 

transformed parameters { 
    real<lower=0 , upper=1> theta[I]; 
    // theta = inv_logit(gam); // Alternatively can use the in-built inv_logit function which is vectorised 
    for(j in 1:I){ 
     theta[j] = 1/(1 + exp(-gam[j])); 
    } 
} 

model { 
    target += gamma_lpdf(lambda | a, b); 
    target += normal_lpdf(mu | m , 1/sqrt(p)); 
    target += normal_lpdf(gam | mu, 1/sqrt(lambda)); 
    target += binomial_lpmf(x | n , theta); 
} 

Answer 1

В качестве подсказки, попробуйте положить gam (ма) в parameters блоке, а затем объявить и определить theta в transformed parameters блоке в соответствии с распределениями, которые Вы дали в самом начале.

Начинающие к Stan часто предполагают, что Стэн наделен способностью логически выработать последствия вашей программы Стан, когда на самом деле он получает transpiled достаточно буквально C++ и строки кода из transformed parameters и model блоков выполняются над и снова.

Причина, по которой имеет значение gam (ma) или theta - это примитивный параметр, связанный с принципом изменения переменных. Если вы действительно этого хотели, вы могли бы получить те же результаты с исходной параметризацией, если бы добавили определяющий термин якобиана (в единицах журнала) до target, но этого легче избежать, переместив gam (ma) в блок parameters и theta к блоку transformed parameters. Подробнее о принципе изменения переменных см. Здесь case study.