将参数固定到 JAGS 中的分布
Fixing a parameter to a distribution in JAGS
在贝叶斯编程语言 JAGS 中，我正在寻找一种方法来将参数固定为特定分布，而不是常数。下面的段落更明确地提出了这个问题并引用了 JAGS 代码。我也愿意接受使用其他概率编程语言（例如 stan）的答案。
下面的第一个代码块 (model1) 是一个 JAGS 脚本，旨在估计具有不等方差的两组高斯混合模型。我正在寻找一种方法来将其中一个参数（例如 $\mu_2$）固定为特定分布（例如，dnorm(0,0.0001)）。我知道如何将 $\mu_2$ 固定为常量（例如，请参阅代码块 2 中的模型 2），但我找不到将 $\mu_2$ 固定为我先前的信念的方法（例如，请参阅代码块 3 中的模型 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
}
"
我不知道 JAGS，但这里有两个 Stan 版本。一个在所有迭代中取 mu2 的单个样本；第二个为每次迭代采用不同的 mu2 样本。
不管怎样，我没有资格判断这是否真的是个好主意。 （特别是第二个版本，是 Stan 团队有意避免的，原因已描述 here。）但它至少是可能的。
（在这两个例子中，我更改了一些先验分布以使数据更易于处理，但基本思想是相同的。）
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 的新样本
我们无法在参数中生成一个运行dom number，t运行sformed parameters, or model block；同样，这是一个深思熟虑的设计选择。但是我们可以在 t运行sformed 数据块中生成一大堆数字，并为每次迭代获取一个新数字。为此，我们需要一种方法来确定我们在参数块中进行的迭代。我从 this discussion on the Stan forums 末尾开始使用 Louis 的解决方案。首先，在您的工作目录中将以下 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();
}
请注意，该模型实际上为 mu2 生成了比我们需要的多 1 个 运行dom 值。当我尝试准确生成 n_iterations 运行dom 值时，我收到一条错误消息，通知我 Stan 已尝试访问 all_mu2s[1001]。
我觉得这很令人担忧，因为这意味着我不完全理解内部发生的事情——考虑到下面的 R 代码，难道不应该只有 1000 次迭代吗？不过看起来就是差一的错误，而且拟合出来的模型看起来还算合理，所以就没再追究了。
另请注意，此方法获取的是迭代次数，而不是链。我运行只有一个链条；如果您 运行 多个链，则 mu2 的第 i 个值在每个链中都是相同的。同一个 Stan 论坛讨论有区分链的建议，但我没有探索它。
最后，生成假数据并将模型与其拟合。当我们编译模型时，我们需要潜入 iter.hpp 中的函数定义，如 here.
所述
# 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 使用了整个 运行ge 个值，它们遵循指定的分布。
# 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")
在贝叶斯编程语言 JAGS 中，我正在寻找一种方法来将参数固定为特定分布，而不是常数。下面的段落更明确地提出了这个问题并引用了 JAGS 代码。我也愿意接受使用其他概率编程语言（例如 stan）的答案。
下面的第一个代码块 (model1) 是一个 JAGS 脚本，旨在估计具有不等方差的两组高斯混合模型。我正在寻找一种方法来将其中一个参数（例如 $\mu_2$）固定为特定分布（例如，dnorm(0,0.0001)）。我知道如何将 $\mu_2$ 固定为常量（例如，请参阅代码块 2 中的模型 2），但我找不到将 $\mu_2$ 固定为我先前的信念的方法（例如，请参阅代码块 3 中的模型 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
}
"
我不知道 JAGS，但这里有两个 Stan 版本。一个在所有迭代中取 mu2 的单个样本；第二个为每次迭代采用不同的 mu2 样本。
不管怎样，我没有资格判断这是否真的是个好主意。 （特别是第二个版本，是 Stan 团队有意避免的，原因已描述 here。）但它至少是可能的。
（在这两个例子中，我更改了一些先验分布以使数据更易于处理，但基本思想是相同的。）
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 的新样本
我们无法在参数中生成一个运行dom number，t运行sformed parameters, or model block；同样，这是一个深思熟虑的设计选择。但是我们可以在 t运行sformed 数据块中生成一大堆数字，并为每次迭代获取一个新数字。为此，我们需要一种方法来确定我们在参数块中进行的迭代。我从 this discussion on the Stan forums 末尾开始使用 Louis 的解决方案。首先，在您的工作目录中将以下 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();
}
请注意，该模型实际上为 mu2 生成了比我们需要的多 1 个 运行dom 值。当我尝试准确生成 n_iterations 运行dom 值时，我收到一条错误消息，通知我 Stan 已尝试访问 all_mu2s[1001]。
我觉得这很令人担忧，因为这意味着我不完全理解内部发生的事情——考虑到下面的 R 代码，难道不应该只有 1000 次迭代吗？不过看起来就是差一的错误，而且拟合出来的模型看起来还算合理，所以就没再追究了。
另请注意，此方法获取的是迭代次数，而不是链。我运行只有一个链条；如果您 运行 多个链，则 mu2 的第 i 个值在每个链中都是相同的。同一个 Stan 论坛讨论有区分链的建议，但我没有探索它。
最后，生成假数据并将模型与其拟合。当我们编译模型时，我们需要潜入 iter.hpp 中的函数定义，如 here.
# 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 使用了整个 运行ge 个值，它们遵循指定的分布。
# 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")