Как решать, когда один параметр зависит от другого в RSTAN - PullRequest
Новая
       66

Как решать, когда один параметр зависит от другого в RSTAN

0 голосов
/ 25 июня 2018

У меня есть код Стэна, в котором один из параметров модели зависит от другого параметра. Всего у меня 5 параметров: му, альфа, бета, гамма, дельта . Теперь бета зависит от альфы таким образом, что - 

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;
}

Как поместить нижнюю границу бета в блок параметров?

Код: 

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)

1 Ответ

0 голосов
/ 25 июня 2018

tldr;

Вы можете реализовать ограничения параметров, основанные на других параметрах, через границы lower / upper.

Простой пример

Давайте создадим простой пример, подгоняя нормальное распределение к некоторым случайным данным, чтобы оценить параметры нормали mu (среднее) и 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 в блоке 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;
...
...