Смешанная модель Pareto и Normal Stan не работает

Question

Я пытаюсь изучить Stan через rstan (так как я знаком с R). Я попытался запустить простую смешанную модель Pareto и Normal. Он компилируется отлично (насколько я могу судить), но он не может пробовать, давая мне ошибку:

«Инициализация между (-2, 2) не удалась после 100 попыток. Попробуйте указать начальные значения, уменьшив диапазоны Зависимые значения, или reparameterizing модели не произошла

ошибки при вызове пробоотборника;. выборки не сделал»

Достаточно сказать, что я пробовал различные способы спараметрировать вещи, и попыталась установить начальные значения, но все нет помогло.

Мой R + rstan код ниже:

library(rstan) 
rpareto = function(n, location, shape){location/runif(n)^(1/shape)} 
sdvec=runif(1e3,0.1,1) 
HMFtest=list(x=rpareto(1e3,10,2)+rnorm(1e3,0,sdvec), sdev=sdvec, N=1e3) 

HMF.stan <- " 
data { 
    int<lower=0> N; 
    real x[N]; 
    real sdev[N]; 
} 
parameters { 
    real<lower=0,upper=20> y_min; 
    real<lower=0,upper=4> alpha; 
    real xtrue[N]; 
} 
model { 
    y_min ~ lognormal(1, 1); 
    alpha ~ lognormal(1, 1); 
    xtrue ~ pareto(y_min, alpha); 
    for(i in 1:N){ 
    x[i] ~ normal(xtrue[i], sdev[i]); 
    } 
} 
" 

stan.test <- stan(model_code=HMF.stan, data=HMFtest, pars=c('y_min','alpha'), chains=1, iter=30000, warmup=10000) 

Этот пример прекрасно работает с зазубринами (следовательно, я меченый Зубцы тоже), и я могу разместить этот код полезно.

Кстати, если я изменю распределение Парето на дополнительное нормальное распределение, он отлично работает (но, конечно, дает мне бессмысленный ответ).

Любые предложения относительно того, что я делаю неправильно, будут высоко оценены! Я боюсь, как-то я все еще думаю, что JAGS не Стэн, но я не мог найти любые примеры людей, подгоняющих модели Pareto со Стэном, поэтому мне было трудно перекреститься, чтобы подтвердить мой подход.

Answer 1

Сообщение об ошибке означает, что все случайные стартовые точки дали вероятность нуля.

Я был в состоянии воспроизвести вашу проблему в Стан, с моделью

data { 
    int<lower=0> N; 
    real x[N]; 
    real sdev[N]; 
} 
parameters { 
    real<lower=0,upper=20> y_min; 
    real<lower=0,upper=4> alpha; 
    real xtrue[N]; 
} 
model { 
    y_min ~ lognormal(1, 1); 
    alpha ~ lognormal(1, 1); 
    print("y_min=", y_min, " alpha=", alpha); 
    xtrue ~ pareto(y_min, alpha); 
    print("xtrue: ", xtrue); 
    x ~ normal(xtrue, sdev); 
    print("x=", x); 
} 

и данные

N <- 6 
sdev <- c(0.3339302,0.2936877,0.8540434,0.2399283,0.1014759,0.3717446) 
x <- c(12.640112,10.502748,11.015629,29.382395,61.180509,12.772482) 

Компиляция и запуск с Stan 2.0.1 (довольно старый сейчас) я получаю выход вроде следующие:

y_min=4.49609:0 alpha=2.54906:0 
xtrue: [0.992331:0,0.303142:0,0.180334:0,1.96009:0,0.903113:0,1.75711:0] 
x=[12.6401,10.5027,11.0156,29.3824,61.1805,12.7725] 
y_min=17.0143:0 alpha=1.67509:0 
xtrue: [-1.40618:0,1.82026:0,1.67344:0,-0.973618:0,0.746502:0,1.93469:0] 
x=[12.6401,10.5027,11.0156,29.3824,61.1805,12.7725] 

Так что, когда разумные параметры выбраны для y_min и alpha, парето- генерируемые значения также ниже y_min. В руководстве функция распределения вероятности также не содержит усечения. Я думаю, что это проблема (замена парето нормальным распределением работает нормально). Я рекомендую открыть ошибку со Стэном на github, заявив, что x ~ pareto (y_min, alpha) генерирует значения ниже y_min.

Код работает с последней версией Stan. Сначала попробуйте обновить, ошибка, похоже, была исправлена ​​некоторое время назад.

Answer 2

Основной проблемой является несоответствие между опорой заявленной на параметры parameters { real<lower=0,upper=20> y_min; real<lower=0,upper=4> alpha; real xtrue[N]; } и образец пространстве априориями model { y_min ~ lognormal(1, 1); alpha ~ lognormal(1, 1); xtrue ~ pareto(y_min, alpha); ...

1. y_min ограничена в (0,20) интервала, но логнормальное предварительно сбрасывает единицу массы по всей положительной реальной линии

2. alpha ограничивается интервалом (0,20), но с ошибкой предшествует разбрасывание единицы массы по всей положительной вещественной линии

3. Худший из всех, каждый элемент xtrue не имеет ограничений - это означает, что ему разрешено быть чем угодно по всей реальной линии --- но ранее парето разбрасывает единицу массы на интервале (y_min, бесконечность)

Проще всего сделать было бы объявить параметры, как parameters { real<lower=0> y_min; real<lower=0> alpha; real<lower=y_min> xtrue[N]; } в принципе, вы можете держать верхние границы y_min и alpha и указать некоторые априорные, которые интегрирорваться 1 над заявленной поддержкой. Грубый способ сделать это состоит в том, чтобы усечь (который будет делить логарифмический PDF на количество необрезанной массы), логнормальные приоритеты, такие как model { y_min ~ lognormal(1, 1) T[,20]; alpha ~ lognormal(1, 1) T[,4]; Возможно, равномерное или четырехпараметрическое бета-распределение было бы более уместным, чем усеченный логарифмический.

Наконец, хотя это не является логически неправильным for(i in 1:N){ x[i] ~ normal(xtrue[i], sdev[i]); } гораздо хуже, чем в вычислительном логически эквивалентное заявление x ~ normal(xtrue, sdev);

Answer 3

Все фиксированные люди. Оказывается, были проблемы глубокие с моей смесью R, stan и компилятора C++.

Не желая беспокоиться об исправлении вещей вручную, я взял подход к кувалде и просто обновил до Йосемити, а затем установил ключевые компоненты с нуля. Похоже, что все это исправляет, и теперь мои цепи сходятся очень красиво.

Это было странно, поскольку компилятор мог построить rstan и скомпилировать/сэмпл многих примеров Rast. Я не знаю, почему парето вызвало эти проблемы, но они, безусловно, исправлены.