当一个参数依赖于RSTAN中的另一个时如何解决

Question

我有一个 Stan 代码，其中一个模型参数取决于另一个参数。我总共有 5 个参数：mu、alpha、beta、gamma、delta。现在 beta 在某种程度上依赖于 alpha-

beta> 1- (alpha/1.17)

.参数块目前的样子：

parameters{
real<lower=0> mu;
real<lower=0,upper=1.17> alpha;
real beta;
real<lower=1> gamma;
real<lower=0> delta;
}

如何将 beta 的下限放入参数块中？

代码是：

expcode="
    functions{
real loglikelihood(int N,
real mu,
real alpha,
real beta,
real gamma,
real delta,
real[] t,
real[] m,
real[] rts,
real magmin,
real tmax,
real betalim){

  real tempA;
  real sumtermA;
  real a;
  real tempB;
  real final;

  sumtermA=log(mu);
  for(j in 2:N){
   tempA=mu;
   for(i in 1:(j-1)){
    tempA += beta*(exp(alpha*(m[i]-magmin)))*(gamma - 1) * delta^(gamma- 1) *(1 / (t[j]-t[I]+delta)^gamma);
   }
   sumtermA += log(tempA);
  }
  tempB=0;
  for(j in 1:N){
   tempB += beta*(exp(alpha*(m[j]-magmin)))*(1-((delta^(gamma- 1))/((tmax-t[j]+delta)^(gamma-1))));
}

 a= mu*tmax;
 final= sumtermA-a-tempB+sum(rts);
 return(final);

  }
}
      data{
        int<lower=0> N;
        real<lower=0> t[N];
        real<lower=0> m[N];
        real rts[N];
        real<lower=0> tmax;
        real<lower=0> magmin;
        real<lower=0> betalim;
        }
       parameters{
     real<lower=0> mu;
     real<lower=0,upper=betalim> alpha;
     real<lower=0> beta;
     real<lower=1> gamma;
     real<lower=0> delta;
     }
     model{
      mu~normal(1.5,1.5);
      alpha~normal(0,0.1);
      beta~normal(0,0.1);
      gamma~normal(1.12,0.16);
      delta~gamma(0.1,0.1);

       //likelihood
          target+=    (loglikelihood(N,mu,alpha,beta,gamma,delta,t,m,rts,magmin,tmax,betalim));
        }
   "

   data<-   list(N=300,t=runif(300,0,1),m=runif(300,2,9),rts=runif(300,-3,3),tmax=1,magmin=2,betalim=1.17)

Answer 1

tldr;

您可以通过lower/upper边界实现基于其他参数的参数约束。

---

一个简单的例子

让我们创建一个简单的示例，通过将正态分布拟合到一些随机数据来估计正态的参数mu（平均值）和sigma（标准差），以及转换后的参数nu = 1/sigma。我们施加约束nu > 1 / sigma - 1。

1. 首先，让我们定义我们的模型。为简单起见，我将使用平面（即默认）先验。

   model <- "
   data {
       int N;                                       // Number of observations
       real y[N];                                   // Response
   }
   
   parameters {
       real mu;                                     // Model parameters
       real<lower=1e-5> sigma;                      // Standard deviation
   }
   
   transformed parameters {
       real<lower = 1 / sigma - 1, upper = positive_infinity()> nu;
       nu = 1 / sigma;
   }
   
   model {
       y ~ normal(mu, sigma);
   }
   "
   

   我们在块transformed parameters 中定义并声明转换后的参数nu。此外，我们通过下限/上限<lower=1/sigma-1, upper=positive_infinity()> 施加约束nu > 1 / sigma - 1。

2. 让我们生成一些示例数据。这里我们选择mu = 2和nu = 1/sigma = 4。

   set.seed(2017);
   mu <- 2;
   nu <- 4;
   y <- rnorm(100, mean = mu, sd = 1/nu);
   

3. 让我们拟合模型。

   library(rstan);
   options(mc.cores = parallel::detectCores())
   rstan_options(auto_write = TRUE)
   fit <- stan(model_code = model, data = list(N = length(y), y = y));
   fit;
   #Inference for Stan model: 16495e5aad9d987998077084a6630917.
   #4 chains, each with iter=2000; warmup=1000; thin=1;
   #post-warmup draws per chain=1000, total post-warmup draws=4000.
   #
   #       mean se_mean   sd  2.5%   25%   50%   75% 97.5% n_eff Rhat
   #mu     2.00    0.00 0.03  1.95  1.98  2.00  2.02  2.06  3382    1
   #sigma  0.27    0.00 0.02  0.24  0.26  0.27  0.28  0.31  3114    1
   #nu     3.72    0.00 0.26  3.21  3.54  3.71  3.89  4.25  3128    1
   #lp__  80.24    0.02 0.96 77.66 79.84 80.54 80.93 81.19  1719    1
   #
   #Samples were drawn using NUTS(diag_e) at Mon Jun 25 11:01:06 2018.
   #For each parameter, n_eff is a crude measure of effective sample size,
   #and Rhat is the potential scale reduction factor on split chains (at
   #convergence, Rhat=1).
   

   您可以看到mu 和nu 的估计值与我们选择的参数值非常一致，而且确实nu > 1 / sigma - 1。

---

你的情况

您应该能够通过将beta 声明为来施加约束beta > 1 - alpha / 1.17

...
real<lower=0,upper=betalim> alpha;
real<lower = 1 - alpha / 1.17, upper = positive_infinity()> beta;
...