В байесовском языке программирования JAGS я ищу способ привязать параметр к определенному распределению, а не к константе. В приведенном ниже параграфе этот вопрос излагается более четко и ссылается на код JAGS. Я также был бы открыт для ответов, в которых используются другие вероятностные языки программирования (например, stan).

Первый фрагмент кода ниже (model1) представляет собой сценарий JAGS, предназначенный для оценки двухгрупповой модели смеси Гаусса с неравными отклонениями. Я ищу способ исправить один из параметров (скажем, $mu_2$) для определенного распределения (например, dnorm(0,0.0001)). Я знаю, как привязать $mu_2$ к константе (например, см. model2 в фрагменте кода 2), хотя я не могу найти способ привязать $mu_2$ к моему предыдущему убеждению(например, см. model3 в фрагменте кода 3, который концептуально показывает, что я пытаюсь сделать).

Заранее спасибо!

Фрагмент кода 1

 model1 = "
model {

for (i in 1:n1){
y1[i] ~ dnorm (mu1 , phi1) 
}

for (i in 1:n2){ 
y2[i] ~ dnorm (mu2 , phi2) 
}

# Priors
phi1 ~ dgamma(.001,.001) 
phi2 ~ dgamma(.001,.001) 
sigma2.1 <- 1/phi1
sigma2.2 <- 1/phi2
mu1 ~ dnorm (0,0.0001) 
mu2 ~ dnorm (0,0.0001)

# Create a variable for the mean difference

delta <- mu1 - mu2

}
"
 

Фрагмент кода 2

 model2 = "
model {

for (i in 1:n1){
y1[i] ~ dnorm (mu1 , phi1) 
}

for (i in 1:n2){ 
y2[i] ~ dnorm (mu2 , phi2) 
}

# Priors
phi1 ~ dgamma(.001,.001) 
phi2 ~ dgamma(.001,.001) 
sigma2.1 <- 1/phi1
sigma2.2 <- 1/phi2
mu1 ~ dnorm (0,0.0001) 
mu2 <- 1.27

# Create a variable for the mean difference

delta <- mu1 - mu2

}
"
 

Фрагмент кода 3

 model3 = "
model {

for (i in 1:n1){
y1[i] ~ dnorm (mu1 , phi1) 
}

for (i in 1:n2){ 
y2[i] ~ dnorm (mu2 , phi2) 
}

# Priors
phi1 ~ dgamma(.001,.001) 
phi2 ~ dgamma(.001,.001) 
sigma2.1 <- 1/phi1
sigma2.2 <- 1/phi2
mu1 ~ dnorm (0,0.0001) 
mu2 <- dnorm (0,0.0001) 

# Create a variable for the mean difference

delta <- mu1 - mu2

}
"
 

Я не знаю ДЖАГОВ, но вот две версии Стэна. Один берет одну выборку mu2 для всех итераций; второй берет другую выборку mu2 для каждой итерации.

В любом случае, я не компетентен судить, действительно ли это хорошая идея. (Вторая версия, в частности, — это то, чего команда Stan намеренно пыталась избежать по причинам, описанным здесь.) Но это, по крайней мере, возможно.

(В обоих примерах я изменил некоторые из предыдущих дистрибутивов, чтобы упростить работу с данными, но основная идея та же.)

Один образец mu2

Во-первых, модель Стэна.

 data {
  int<lower=0> n1;
  vector[n1] y1;
  int<lower=0> n2;
  vector[n2] y2;
}

transformed data {
  // Set mu2 to a single randomly selected value (instead of giving it a prior
  // and estimating it).
  real mu2 = normal_rng(0, 0.0001);
}

parameters {
  real mu1;
  real<lower=0> phi1;
  real<lower=0> phi2;
}

transformed parameters {
  real sigma1 = 1 / phi1;
  real sigma2 = 1 / phi2;
}

model {
  mu1 ~ normal(0, 0.0001);
  phi1 ~ gamma(1, 1);
  phi2 ~ gamma(1, 1);
  y1 ~ normal(mu1, sigma1);
  y2 ~ normal(mu2, sigma2);
}

generated quantities {
  real delta = mu1 - mu2;
  // We can't return mu2 from the transformed data block.  So if we want to see
  // what it was, we have to copy its value into a generated quantity and return
  // that.
  real mu2_return = mu2;
}
 

Затем введите код R для генерации поддельных данных и подгонки под модель.

 # Generate fake data.
n1 = 1000
n2 = 1000
mu1 = rnorm(1, 0, 0.0001)
mu2 = rnorm(1, 0, 0.0001)
phi1 = rgamma(1, shape = 1, rate = 1)
phi2 = rgamma(1, shape = 1, rate = 1)
y1 = rnorm(n1, mu1, 1 / phi1)
y2 = rnorm(n2, mu2, 1 / phi2)
delta = mu1 - mu2

# Fit the Stan model.
library(rstan)
options(mc.cores = parallel::detectCores())
rstan_options(auto_write = T)

stan.data = list(n1 = n1, y1 = y1, n2 = n2, y2 = y2)
stan.model = stan(file = "stan_model.stan",
                  data = stan.data,
                  cores = 3, iter = 1000)
 

Мы можем извлечь образцы из модели Stan и убедиться, что мы правильно восстановили истинные значения параметров — за исключением, конечно, случая mu2 .

 # Pull out the samples.
library(tidybayes)
library(tidyverse)
stan.model %>%
  spread_draws(mu1, phi1, mu2_return, phi2) %>%
  ungroup() %>%
  dplyr::select(.draw, mu1, phi1, mu2 = mu2_return, phi2) %>%
  pivot_longer(cols = -c(.draw), names_to = "parameter") %>%
  ggplot(aes(x = value))  
  geom_histogram()  
  geom_vline(data = data.frame(parameter = c("mu1", "phi1", "mu2", "phi2"),
                               true.value = c(mu1, phi1, mu2, phi2)),
             aes(xintercept = true.value), color = "red", size = 1.5)  
  facet_wrap(~ parameter, scales = "free")  
  theme_bw()  
  scale_x_continuous("Parameter value")  
  scale_y_continuous("Number of samples")
 

Новый образец mu2 для каждой итерации

Мы не можем генерировать случайное число в параметрах, преобразованных параметрах или блоке модели; опять же, это преднамеренный выбор дизайна. Но мы можем сгенерировать целую кучу чисел в преобразованном блоке данных и получать новое для каждой итерации. Для этого нам нужен способ определить, на какой итерации мы находимся в блоке параметров. Я использовал решение Луиса с конца этого обсуждения на форумах Stan. Во-первых, сохраните следующий код C как iter.hpp в вашем рабочем каталоге:

 static int itct = 1;
inline void add_iter(std::ostream* pstream__) {
    itct  = 1;
}
inline int get_iter(std::ostream* pstream__) {
    return itct;
}
 

Затем определите модель Stan следующим образом. Функции add_iter() и get_iter() определены в iter.hpp ; если вы работаете в RStudio, вы получите символы ошибок при редактировании файла Stan, потому что RStudio не знает, что мы собираемся ввести эти определения функций из других источников.

 functions {
  void add_iter();
  int get_iter();
}

data {
  int<lower=0> n1;
  vector[n1] y1;
  int<lower=0> n2;
  vector[n2] y2;
  int<lower=0> n_iterations;
}

transformed data {
  vector[n_iterations   1] all_mu2s;
  for(n in 1:(n_iterations   1)) {
    all_mu2s[n] = normal_rng(0, 0.0001);
  }
}

parameters {
  real mu1;
  real<lower=0> phi1;
  real<lower=0> phi2;
}

transformed parameters {
  real sigma1 = 1 / phi1;
  real sigma2 = 1 / phi2;
  real mu2 = all_mu2s[get_iter()];
}

model {
  mu1 ~ normal(0, 0.0001);
  phi1 ~ gamma(1, 1);
  phi2 ~ gamma(1, 1);
  y1 ~ normal(mu1, sigma1);
  y2 ~ normal(mu2, sigma2);
}

generated quantities {
  real delta = mu1 - mu2;
  add_iter();
}
 

Обратите внимание, что модель фактически генерирует на 1 случайное значение больше mu2 , чем нам нужно. Когда я попытался сгенерировать точно n_iterations случайные значения, я получил сообщение об ошибке, сообщающее мне, что Стэн пытался получить доступ all_mu2s[1001] .
Я нахожу это тревожным, потому что это означает, что я не до конца понимаю, что происходит внутри — разве не должно быть только 1000 итераций, учитывая приведенный ниже код R? Но это просто выглядит как разовая ошибка, и подходящая модель выглядит разумно, поэтому я не стал развивать эту тему дальше.

Кроме того, обратите внимание, что этот подход получает номер итерации, но не цепочку. Я запустил только одну цепочку; если вы запустите более одной цепочки, i-е значение mu2 будет одинаковым в каждой цепочке. В том же обсуждении на форумах Стэна есть предложение о различении цепочек, но я его не изучал.

Наконец, создайте поддельные данные и подгоните к ним модель. Когда мы компилируем модель, нам нужно проникнуть в определения функций iter.hpp , как описано здесь.

 # Generate fake data.
n1 = 1000
n2 = 1000
mu1 = rnorm(1, 0, 0.0001)
mu2 = rnorm(1, 0, 0.0001)
phi1 = rgamma(1, shape = 1, rate = 1)
phi2 = rgamma(1, shape = 1, rate = 1)
y1 = rnorm(n1, mu1, 1 / phi1)
y2 = rnorm(n2, mu2, 1 / phi2)
delta = mu1 - mu2
n.iterations = 1000

# Fit the Stan model.
library(rstan)
stan.data = list(n1 = n1, y1 = y1, n2 = n2, y2 = y2,
                 n_iterations = n.iterations)
stan.model = stan_model(file = "stan_model.stan",
                        allow_undefined = T,
                        includes = paste0('n#include "',
                                          file.path(getwd(), 'iter.hpp'),
                                          '"n'))
stan.model.fit = sampling(stan.model,
                          data = stan.data,
                          chains = 1,
                          iter = n.iterations,
                          pars = c("mu1", "phi1", "mu2", "phi2"))
 

В очередной раз мы восстановили ценности mu1 phi1 , и phi2 достаточно хорошо. На этот раз мы использовали целый диапазон значений для mu2 , которые следуют указанному распределению.

 # Pull out the samples.
library(tidybayes)
library(tidyverse)
stan.model.fit %>%
  spread_draws(mu1, phi1, mu2, phi2) %>%
  ungroup() %>%
  dplyr::select(.draw, mu1, phi1, mu2 = mu2, phi2) %>%
  pivot_longer(cols = -c(.draw), names_to = "parameter") %>%
  ggplot(aes(x = value))  
  geom_histogram()  
  stat_function(dat = data.frame(parameter = "mu2", value = 0),
                fun = function(.x) { dnorm(.x, 0, 0.0001) * 0.01 },
                color = "blue", size = 1.5)  
  geom_vline(data = data.frame(parameter = c("mu1", "phi1", "mu2", "phi2"),
                               true.value = c(mu1, phi1, mu2, phi2)),
             aes(xintercept = true.value), color = "red", size = 1.5)  
  facet_wrap(~ parameter, scales = "free")  
  theme_bw()  
  scale_x_continuous("Parameter value")  
  scale_y_continuous("Number of samples")