需要修复广义帕累托分布的 Stan 代码以接受实参而不是向量
Need to fix Stan code for the generalized pareto distribution to take in real arguments rather than vectors
我正在使用此处定义的函数：Stan 中的极值分析和用户定义的概率函数，用于使用广义帕累托分布对数据进行建模，但我的问题是我的模型处于 for 循环中，需要三个实数值参数，而 gpd 函数假设一个向量、实数、实数参数。
我不太确定我的模型块是否适合矢量化，所以我想我需要让 gpd 函数接受实值参数（但也许我错了）。
如果能帮助我切换代码以实现此目的，我将不胜感激。这是我的斯坦码
functions {
real gpareto_lpdf(vector y, real k, real sigma) {
// generalised Pareto log pdf
int N = rows(y);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y) * (k/sigma))) -N*log(sigma);
else
return -sum(y)/sigma -N*log(sigma); // limit k->0
}
real gpareto_lcdf(vector y, real k, real sigma) {
// generalised Pareto log cdf
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return sum(log1m_exp((-inv_k)*(log1p((y) * (k/sigma)))));
else
return sum(log1m_exp(-(y)/sigma)); // limit k->0
}
}
data {
// the input data
int<lower = 1> n; // number of observations
real<lower = 0> value[n]; // value measurements
int<lower = 0, upper = 1> censored[n]; // vector of 0s and 1s
// parameters for the prior
real<lower = 0> a;
real<lower = 0> b;
}
parameters {
real k;
real sigma;
}
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
for (i in 1:n) {
if (censored[i]) {
target += gpareto_lcdf(value[i] | k, sigma);
} else {
target += gpareto_lpdf(value[i] | k, sigma);
}
}
}
这是日志 PDF 的调整方式。这样，可以传递用于将 y 子集化为截尾和 non-censored 观察的索引数组。
real cens_gpareto_lpdf(vector y, int[] cens, int[] no_cens, real k, real sigma) {
// generalised Pareto log pdf
int N = size(cens);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y[no_cens]) * (k/sigma))) -N*log(sigma) +
sum(log1m_exp((-inv_k)*(log1p((y[cens]) * (k/sigma)))));
else
return -sum(y[no_cens])/sigma -N*log(sigma) +
sum(log1m_exp(-(y[cens])/sigma));
}
扩展 data 块：n_cens、n_not_cens、cens 和 no_cens 是需要提供的值。
int<lower = 1> n; // total number of obs
int<lower = 1> n_cens; // number of censored obs
int<lower = 1> n_not_cens; // number of regular obs
int cens[n_cens]; // index set censored
int no_cens[n_not_cens]; // index set regular
vector<lower = 0>[n] value; // value measurements
gfgm 建议的非零参数：
parameters {
real<lower=0> k;
real<lower=0> sigma;
}
重写 model 块：
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
value ~ cens_gpareto(cens, no_cens, k, sigma);
}
免责声明：我既没有检查公式的完整性，也没有运行使用测试数据检查模型。刚刚通过 rstan::stan_model() 编译，效果很好。 gfgm 的建议对于 post-processing / generated quantities 等中的计算内容可能更方便。我不是 Stan 专家:-)。
编辑：
修复了 gfgm 通过模拟发现的分歧问题。可能性是 ill-defined（N= rows(y) 而不是 N=size(cens)。现在使用 gfgm 的数据运行良好（使用 set.seed(123) 和 rstan）：
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
k 0.16 0.00 0.10 0.02 0.08 0.14 0.21 0.42 1687 1
sigma 0.90 0.00 0.12 0.67 0.82 0.90 0.99 1.16 1638 1
lp__ -106.15 0.03 1.08 -109.09 -106.56 -105.83 -105.38 -105.09 1343 1
您可以让模型接受向量并避免 for 循环。在某种程度上，您需要在声明 value 时更改签名，但随后还将经过审查和未经审查的观察结果的索引作为数据输入。
我在下方发布 运行 合成数据的代码，但我合成的数据只是一些随机 log-normal 变量，而不是帕累托的实际绘制，所以我没有了解模型是否正在恢复参数或它的覆盖范围是什么样的。您可能需要进行一些基于模拟的校准来检查模型。
斯坦码：
functions {
real gpareto_lpdf(vector y, real k, real sigma) {
// generalised Pareto log pdf
int N = rows(y);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y) * (k/sigma))) -N*log(sigma);
else
return -sum(y)/sigma -N*log(sigma); // limit k->0
}
real gpareto_lcdf(vector y, real k, real sigma) {
// generalised Pareto log cdf
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(k) > 1e-15)
return sum(log1m_exp((-inv_k)*(log1p((y) * (k/sigma)))));
else
return sum(log1m_exp(-(y)/sigma)); // limit k->0
}
}
data {
// the input data
int<lower = 1> n; // number of observations
int<lower = 0> n_cens;
vector<lower = 0>[n] value; // value measurements
int<lower=1> cens_id[n_cens]; // the indices of censored observations
int<lower=1> nocens_id[n - n_cens]; // the indices of uncensored obs
// parameters for the prior
real<lower = 0> a;
real<lower = 0> b;
}
parameters {
real<lower=0> k;
real<lower=0> sigma;
}
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
target += gpareto_lcdf(value[cens_id] | k, sigma);
target += gpareto_lpdf(value[nocens_id] | k, sigma);
}
这 运行 在我的 mid-tier 笔记本电脑上使用 cmdstanr 基本上是即时的。您可以 运行 使用 R 代码
library(cmdstanr)
library(bayesplot)
gparet_mod <- cmdstan_model("gpareto_example.stan")
# make some fake data. This won't be correct and its
# no proper test of the capacity of the model to
# recover the parameters, just seeing if it runs
N <- 100
value <- rlnorm(N)
censored <- rbinom(N, 1, .5)
N_cens <- sum(censored)
cens_id <- which(censored == 1)
nocens_id <- which(censored == 0)
a <- b <- 2
dat <- list(n = N, n_cens = N_cens, value = value,
cens_id = cens_id, nocens_id = nocens_id,
a = a, b = b)
ests <- gparet_mod$sample(data = dat, parallel_chains = 2, chains = 2)
这会产生看似正常的样本
ests$summary()
# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -103. -103. 1.10 0.767 -105. -102. 1.00 750. 807.
2 k 0.316 0.296 0.135 0.131 0.128 0.557 1.00 665. 570.
3 sigma 0.691 0.685 0.113 0.110 0.511 0.887 1.00 772. 887.
和后验分布
mcmc_hist(ests$draws(variables = c("k", "sigma")))
我认为 Martin Arnold 的方法更优雅一些，但是当我试图在相同的数据上将其变为 运行 时，采样器因分歧而崩溃，所以也许我移植了错误的模型或者有些东西需要调整。
编辑
鉴于似然函数，将参数 k 和 sigma 限制为正值更有意义。该模型 运行 更好，您会收到更少的错误消息。我修改了上面的 stan 代码以反映这一点。
我正在使用此处定义的函数：Stan 中的极值分析和用户定义的概率函数，用于使用广义帕累托分布对数据进行建模，但我的问题是我的模型处于 for 循环中，需要三个实数值参数，而 gpd 函数假设一个向量、实数、实数参数。
我不太确定我的模型块是否适合矢量化，所以我想我需要让 gpd 函数接受实值参数（但也许我错了）。
如果能帮助我切换代码以实现此目的，我将不胜感激。这是我的斯坦码
functions {
real gpareto_lpdf(vector y, real k, real sigma) {
// generalised Pareto log pdf
int N = rows(y);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y) * (k/sigma))) -N*log(sigma);
else
return -sum(y)/sigma -N*log(sigma); // limit k->0
}
real gpareto_lcdf(vector y, real k, real sigma) {
// generalised Pareto log cdf
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma)
if (fabs(k) > 1e-15)
return sum(log1m_exp((-inv_k)*(log1p((y) * (k/sigma)))));
else
return sum(log1m_exp(-(y)/sigma)); // limit k->0
}
}
data {
// the input data
int<lower = 1> n; // number of observations
real<lower = 0> value[n]; // value measurements
int<lower = 0, upper = 1> censored[n]; // vector of 0s and 1s
// parameters for the prior
real<lower = 0> a;
real<lower = 0> b;
}
parameters {
real k;
real sigma;
}
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
for (i in 1:n) {
if (censored[i]) {
target += gpareto_lcdf(value[i] | k, sigma);
} else {
target += gpareto_lpdf(value[i] | k, sigma);
}
}
}
这是日志 PDF 的调整方式。这样，可以传递用于将 y 子集化为截尾和 non-censored 观察的索引数组。
real cens_gpareto_lpdf(vector y, int[] cens, int[] no_cens, real k, real sigma) {
// generalised Pareto log pdf
int N = size(cens);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma)
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y[no_cens]) * (k/sigma))) -N*log(sigma) +
sum(log1m_exp((-inv_k)*(log1p((y[cens]) * (k/sigma)))));
else
return -sum(y[no_cens])/sigma -N*log(sigma) +
sum(log1m_exp(-(y[cens])/sigma));
}
扩展 data 块：n_cens、n_not_cens、cens 和 no_cens 是需要提供的值。
int<lower = 1> n; // total number of obs
int<lower = 1> n_cens; // number of censored obs
int<lower = 1> n_not_cens; // number of regular obs
int cens[n_cens]; // index set censored
int no_cens[n_not_cens]; // index set regular
vector<lower = 0>[n] value; // value measurements
gfgm 建议的非零参数：
parameters {
real<lower=0> k;
real<lower=0> sigma;
}
重写 model 块：
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
value ~ cens_gpareto(cens, no_cens, k, sigma);
}
免责声明：我既没有检查公式的完整性，也没有运行使用测试数据检查模型。刚刚通过 rstan::stan_model() 编译，效果很好。 gfgm 的建议对于 post-processing / generated quantities 等中的计算内容可能更方便。我不是 Stan 专家:-)。
编辑：
修复了 gfgm 通过模拟发现的分歧问题。可能性是 ill-defined（N= rows(y) 而不是 N=size(cens)。现在使用 gfgm 的数据运行良好（使用 set.seed(123) 和 rstan）：
mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
k 0.16 0.00 0.10 0.02 0.08 0.14 0.21 0.42 1687 1
sigma 0.90 0.00 0.12 0.67 0.82 0.90 0.99 1.16 1638 1
lp__ -106.15 0.03 1.08 -109.09 -106.56 -105.83 -105.38 -105.09 1343 1
您可以让模型接受向量并避免 for 循环。在某种程度上，您需要在声明 value 时更改签名，但随后还将经过审查和未经审查的观察结果的索引作为数据输入。
我在下方发布 运行 合成数据的代码，但我合成的数据只是一些随机 log-normal 变量，而不是帕累托的实际绘制，所以我没有了解模型是否正在恢复参数或它的覆盖范围是什么样的。您可能需要进行一些基于模拟的校准来检查模型。
斯坦码：
functions {
real gpareto_lpdf(vector y, real k, real sigma) {
// generalised Pareto log pdf
int N = rows(y);
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(k) > 1e-15)
return -(1+inv_k)*sum(log1p((y) * (k/sigma))) -N*log(sigma);
else
return -sum(y)/sigma -N*log(sigma); // limit k->0
}
real gpareto_lcdf(vector y, real k, real sigma) {
// generalised Pareto log cdf
real inv_k = inv(k);
if (k<0 && max(y)/sigma > -inv_k)
reject("k<0 and max(y)/sigma > -1/k; found k, sigma =", k, sigma);
if (sigma<=0)
reject("sigma<=0; found sigma =", sigma);
if (fabs(k) > 1e-15)
return sum(log1m_exp((-inv_k)*(log1p((y) * (k/sigma)))));
else
return sum(log1m_exp(-(y)/sigma)); // limit k->0
}
}
data {
// the input data
int<lower = 1> n; // number of observations
int<lower = 0> n_cens;
vector<lower = 0>[n] value; // value measurements
int<lower=1> cens_id[n_cens]; // the indices of censored observations
int<lower=1> nocens_id[n - n_cens]; // the indices of uncensored obs
// parameters for the prior
real<lower = 0> a;
real<lower = 0> b;
}
parameters {
real<lower=0> k;
real<lower=0> sigma;
}
model {
// prior
k ~ gamma(a, b);
sigma ~ gamma(a,b);
// likelihood
target += gpareto_lcdf(value[cens_id] | k, sigma);
target += gpareto_lpdf(value[nocens_id] | k, sigma);
}
这 运行 在我的 mid-tier 笔记本电脑上使用 cmdstanr 基本上是即时的。您可以 运行 使用 R 代码
library(cmdstanr)
library(bayesplot)
gparet_mod <- cmdstan_model("gpareto_example.stan")
# make some fake data. This won't be correct and its
# no proper test of the capacity of the model to
# recover the parameters, just seeing if it runs
N <- 100
value <- rlnorm(N)
censored <- rbinom(N, 1, .5)
N_cens <- sum(censored)
cens_id <- which(censored == 1)
nocens_id <- which(censored == 0)
a <- b <- 2
dat <- list(n = N, n_cens = N_cens, value = value,
cens_id = cens_id, nocens_id = nocens_id,
a = a, b = b)
ests <- gparet_mod$sample(data = dat, parallel_chains = 2, chains = 2)
这会产生看似正常的样本
ests$summary()
# A tibble: 3 × 10
variable mean median sd mad q5 q95 rhat ess_bulk ess_tail
<chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
1 lp__ -103. -103. 1.10 0.767 -105. -102. 1.00 750. 807.
2 k 0.316 0.296 0.135 0.131 0.128 0.557 1.00 665. 570.
3 sigma 0.691 0.685 0.113 0.110 0.511 0.887 1.00 772. 887.
和后验分布
mcmc_hist(ests$draws(variables = c("k", "sigma")))
我认为 Martin Arnold 的方法更优雅一些，但是当我试图在相同的数据上将其变为 运行 时，采样器因分歧而崩溃，所以也许我移植了错误的模型或者有些东西需要调整。
编辑
鉴于似然函数，将参数 k 和 sigma 限制为正值更有意义。该模型 运行 更好，您会收到更少的错误消息。我修改了上面的 stan 代码以反映这一点。