Стандартная версия модели JAGS, которая включает сумму дискретных значений - возможно ли это?

Я пытался запустить эту модель в Стэне. У меня есть работающая версия JAGS (которая возвращает сильно автокоррелированные параметры), и я знаю, как сформулировать ее как CDF двойной экспоненты (с двумя скоростями), которая, вероятно, будет работать без проблем. Однако я хотел бы использовать эту версию в качестве отправной точки для аналогичных, но более сложных моделей.

К настоящему времени у меня есть подозрение, что такая модель в Стэне невозможна. Возможно, из-за дискретности, которую вводит сумма логического значения, Стэн не сможет вычислять градиенты.

Кто-нибудь знает, так ли это, или я что-то еще не так делаю в этой модели? Я вставляю ошибки, которые получаю, под кодом модели.

Большое спасибо заранее Ян

Model:

data {
    int y[11]; 
    int reps[11];
    real soas[11]; 

}
parameters {
    real<lower=0.001,upper=0.200> v1;
    real<lower=0.001,upper=0.200> v2;

}


model {
    real dif[11,96];
    real cf[11];

    real p[11];

    real t1[11,96];
    real t2[11,96];

    for (i in 1:11){
        for (r in 1:reps[i]){     
            t1[i,r]  ~ exponential(v1);
            t2[i,r]  ~ exponential(v2);
            dif[i,r] <-  (t1[i,r]+soas[i]<=(t2[i,r]));

            }
        cf[i] <- sum(dif[i]);
        p[i]  <-cf[i]/reps[i];
        y[i] ~ binomial(reps[i],p[i]); 
    }

}

Вот некоторые фиктивные данные:

psy_dat = { 
         'soas' :  numpy.array(range(-100,101,20)),
            'y' :  [47, 46, 62, 50, 59, 47, 36, 13, 7, 2, 1],
         'reps' :  [48, 48, 64, 64, 92, 92, 92, 64, 64, 48, 48]
          }

А вот ошибки:

DIAGNOSTIC(S) FROM PARSER:
Warning (non-fatal): Left-hand side of sampling statement (~) contains a non-linear transform of a parameter or local variable.
You must call increment_log_prob() with the log absolute determinant of the Jacobian of the transform.
Sampling Statement left-hand-side expression:
get_base1(get_base1(t1,i,"t1",1),r,"t1",2) ~ exponential_log(...)
Warning (non-fatal): Left-hand side of sampling statement (~) contains a non-linear transform of a parameter or local variable.
You must call increment_log_prob() with the log absolute determinant of the Jacobian of the transform.
Sampling Statement left-hand-side expression:
get_base1(get_base1(t2,i,"t2",1),r,"t2",2) ~ exponential_log(...)

И во время выполнения:

Informational Message: The current Metropolis proposal is about to be rejected because of the following issue:
 stan::prob::exponential_log(N4stan5agrad3varE): Random variable is nan:0, but must not be nan!
 If this warning occurs sporadically, such as for highly constrained variable types like covariance matrices, then the sampler is fine,
but if this warning occurs often then your model may be either severely ill-conditioned or misspecified.
 Rejecting proposed initial value with zero density.


Initialization between (-2, 2) failed after 100 attempts. 
 Try specifying initial values, reducing ranges of constrained values, or   reparameterizing the model

Вот рабочая версия этой модели JAGS:

   model {
   for ( n in 1 : N  ) { 
     for (r in 1 : reps[n]){
       t1[r,n] ~ dexp(v1)
       t2[r,n] ~ dexp(v2)
       c[r,n] <- (1.0*((t1[r,n]+durs[n])<=t2[r,n]))
     } 
     p[n] <- max((min(sum(c[,n]) /  (reps[n]),0.99999999999999)),   1-0.99999999999999)) 
     y[n] ~ dbin(p[n],reps[n])
   }

   v1 ~ dunif(0.0001,0.2)
   v2 ~ dunif(0.0001,0.2)
   }

Что касается min () и max (): см. Этот пост https://stats.stackexchange.com/questions/130978/observed-node-inconsistent-when-binomial-success-rate-exactly-one?noredirect=1#comment250046_130978 .


person Jan Tünnermann    schedule 15.01.2015    source источник


Ответы (1)


Я все еще не уверен, какую модель вы пытаетесь оценить (было бы лучше, если вы разместите код JAGS), но то, что у вас есть выше, не может работать в Stan. Стэн ближе к C ++ в том смысле, что вам нужно объявлять, а затем определять объекты. В вашей программе Stan у вас есть два объявления real t1[11,96]; real t2[11,96]; , но нет определений t1 или t2. Следовательно, они инициализируются в NaN, и когда вы выполняете t1[i,r] ~ exponential(v1); , он анализируется как что-то вроде for(i in 1:11) for(r in 1:reps[i]) lp__ += log(v1) - v1 * NaN , где lp__ - это внутренний символ, который содержит значение log-posterior, которое становится NaN и не может обновлять параметры в стиле Metropolis.

Возможно, вы имели в виду, что t1 и t2 должны быть неизвестными параметрами, и в этом случае они должны быть объявлены в блоке parameters. Следующая [EDITED] программа Stan действительна и должна работать, но это может быть не та программа, которую вы имели в виду (для меня это не имеет большого смысла, и прерывность в dif, вероятно, помешает Стэну выполнить выборку из этого апостериорного распределения эффективно). data { int<lower=1> N; int y[N]; int reps[N]; real soas[N]; } parameters { real<lower=0.001,upper=0.200> v1; real<lower=0.001,upper=0.200> v2; real t1[N,max(reps)]; real t2[N,max(reps)]; } model { for (i in 1:N) { real dif[reps[i]]; for (r in 1:reps[i]) { dif[r] <- (t1[i,r]+soas[i]) <= t2[i,r]; } y[i] ~ binomial(reps[i], (1.0 + sum(dif)) / (1.0 + reps[i])); } to_array_1d(t1) ~ exponential(v1); to_array_1d(t2) ~ exponential(v2); }

person Ben Goodrich    schedule 16.01.2015
comment
Я пробовал именно это (и проверял, используя ваши точные предложения): это дает ту же ошибку. Я добавляю код JAGS к исходному вопросу. - person Jan Tünnermann; 16.01.2015
comment
Я добавляю код JAGS к исходному вопросу. Обратите внимание, что я также пробовал использовать функцию min () / max () в Stan. Что касается моделирования: это вероятность того, что один процесс завершится раньше, чем другой. Что, если выразить аналитически, очень хорошо оценивается в JAGS, но в этой форме цепи сильно автокоррелированы. Я надеялся, что HMC Стэна улучшит это. Причина такого выражения модели заключается в том, что она близка к физически предполагаемым процессам и позволит очень легко включить некоторые более продвинутые механизмы без получения сложного аналитического решения. - person Jan Tünnermann; 16.01.2015
comment
В ПОРЯДКЕ. Откуда 96 в исходном коде? Возможно, вы имели в виду 92, это максимальное количество повторений (я изменил код Стэна, чтобы отразить это предположение). Раньше, когда длина dif[i] составляла 96 элементов, и вы выполняли sum(dif[i]), когда reps[i] < 96, вы получали NaN, потому что последние (несколько) элементов dif[i] остаются неинициализированными. Однако это все равно не сработает со Стэном из-за разногласий. Разве биномиальная вероятность не может быть выражена аналитически как CDF сдвинутой экспоненты с верхней границей t2[i,r]? - person Ben Goodrich; 16.01.2015
comment
96 должно быть максимумом (повторений) - глупая ошибка из-за жестко запрограммированных значений. И это был хороший намек. После исправления ошибка осталась, но я заметил, что сумма (почти) всегда равна nan, потому что diff содержит nans всякий раз, когда reps [i] ‹max (reps), что имеет место для многих i. Для тестирования я использовал 92 для повторений [i] при всех i, и сейчас действительно проводится выборка модели. И да, это может быть успешно выражено (по крайней мере, в JAGS) с использованием функций CDF (например, CDF с двойным выражением для моделирования обоих процессов); но посмотрите исходный пост о моей мотивации для этой версии. Вернусь после того, как посмотрю на задние части ... - person Jan Tünnermann; 17.01.2015