0

我有一个 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)
4

1 回答 1

4

tldr;

lower通过/ upperbounds可以实现基于其他参数的参数约束。

一个简单的例子

让我们通过将正态分布拟合到一些随机数据来创建一个简单的示例，以估计正态的参数mu(mean) 和sigma(stdev)，以及转换后的参数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);
}
"

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

  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 > 1 - alpha / 1.17通过声明beta为来施加约束

...
real<lower=0,upper=betalim> alpha;
real<lower = 1 - alpha / 1.17, upper = positive_infinity()> beta;
...
于 2018-06-25T01:12:53.287 回答