在 rstan 中转换变量(贝叶斯分析)
Transforming variables in rstan (bayesian analysis)
我是贝叶斯分析的新手,正在尝试使用 rstan 来估计后验密度分布。该练习试图使用 stan 重新创建我的大学给我们的示例,但我对如何正确转换变量有点困惑。我当前的代码运行没有错误,但结果与大学给我们的结果不符(尽管很接近),为清楚起见,下图使用黑色的 stan 估计。我通过查阅手册并将随机位拼凑在一起来使代码工作,但特别是我不太确定为什么需要 target
以及 gamma 转换是否确实正确。任何指导将不胜感激!
型号
斯坦码
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
)
dat
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")
)
---更新解决方案---
Ben Goodrich 指出的问题是,我是从 theta 推导伽玛的,因为伽玛是我的随机变量,所以它应该是相反的。正确的 stan 代码如下。
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 );
}
作为提示,尝试将 gam
(ma) 放在 parameters
块中,然后根据您提供的分布在 transformed parameters
块中声明和定义 theta
一开始。
Stan 的初学者通常认为 Stan 有能力从逻辑上计算出你的 Stan 程序的含义,但实际上它被相当字面地转译为 C++ 以及来自 transformed parameters
和 model
块被一遍又一遍地执行。
之所以gam
(ma)或theta
是原始参数有区别,与change-of-variables原理有关。如果您真的想要,如果您将雅可比行列式项(以对数单位)添加到 target
,您可以获得与原始参数化相同的结果,但是通过移动 gam
( ma) 到 parameters
块,theta
到 transformed parameters
块。 change-of-variables原理详见case study.
我是贝叶斯分析的新手,正在尝试使用 rstan 来估计后验密度分布。该练习试图使用 stan 重新创建我的大学给我们的示例,但我对如何正确转换变量有点困惑。我当前的代码运行没有错误,但结果与大学给我们的结果不符(尽管很接近),为清楚起见,下图使用黑色的 stan 估计。我通过查阅手册并将随机位拼凑在一起来使代码工作,但特别是我不太确定为什么需要 target
以及 gamma 转换是否确实正确。任何指导将不胜感激!
型号
斯坦码
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
)
dat
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")
)
---更新解决方案---
Ben Goodrich 指出的问题是,我是从 theta 推导伽玛的,因为伽玛是我的随机变量,所以它应该是相反的。正确的 stan 代码如下。
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 );
}
作为提示,尝试将 gam
(ma) 放在 parameters
块中,然后根据您提供的分布在 transformed parameters
块中声明和定义 theta
一开始。
Stan 的初学者通常认为 Stan 有能力从逻辑上计算出你的 Stan 程序的含义,但实际上它被相当字面地转译为 C++ 以及来自 transformed parameters
和 model
块被一遍又一遍地执行。
之所以gam
(ma)或theta
是原始参数有区别,与change-of-variables原理有关。如果您真的想要,如果您将雅可比行列式项(以对数单位)添加到 target
,您可以获得与原始参数化相同的结果,但是通过移动 gam
( ma) 到 parameters
块,theta
到 transformed parameters
块。 change-of-variables原理详见case study.