Гиперприоры для иерархических моделей со Стэном

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

Я видел этот пример для определения гиперприора в pymc, но я не уверен, как сделать что-то подобное со Стэном.

@pymc.stochastic(dtype=np.float64)
def beta_priors(value=[1.0, 1.0]):
    a, b = value
    if a <= 0 or b <= 0:
        return -np.inf
    else:
        return np.log(np.power((a + b), -2.5))

a = beta_priors[0]
b = beta_priors[1]

Затем a и b используются в качестве параметров для предыдущей бета-версии.

Может ли кто-нибудь подсказать мне, как нечто подобное можно было бы сделать со Стэном?


bayesian stan pymc
person Ian_Fin    schedule 12.09.2017    source источник
comment
В руководстве много примеров построения иерархических и многоуровневых моделей (о чем вы здесь говорите). Наши основные рекомендации для априорных задач приведены в главе руководства по регрессии, а также на этой вики-странице: github.com/stan-dev/stan/wiki/Prior-Choice-Recommendations   -  person Bob Carpenter    schedule 13.09.2017
comment
Существуют также тематические исследования по иерархическим моделям, в частности, непосредственно о бинарных переменных, которые противопоставляют гиперприоры для биномов с логистической регрессией только с перехватом (в результате вы, вероятно, не хотите использовать бета-биномы или полиномы Дирихле) : mc-stan.org/users/documentation/case -studies /   -  person Bob Carpenter    schedule 13.09.2017
comment
@BobCarpenter Это действительно полезно - спасибо. Я делаю это для AB-тестирования, а для более сложных, я действительно думал, что работа в регрессионной структуре имеет больше смысла. Возможность ранжировать варианты, как в тематическом исследовании, также может быть полезна для передачи результатов коллегам.   -  person Ian_Fin    schedule 13.09.2017

Ответы (2)


arrow_upward
3
arrow_downward

Чтобы правильно нормализовать это, вам понадобится распределение Парето. Например, если вам нужен дистрибутив p(a, b) ∝ (a + b)^(-2.5), вы можете использовать

a + b ~ pareto(L, 1.5);

где a + b > L. Невозможно нормализовать плотность с поддержкой всех значений, больших или равных нулю - для этого требуется конечное L в качестве нижней границы. Обсуждается использование только этого априорного значения в качестве компонента подсчета иерархического априорного значения для симплекса.

Если a и b являются параметрами, их можно либо ограничить как положительные, либо вы можете оставить a без ограничений и объявить

real<lower = L - a> b;

застраховать a + b > L. L может быть небольшой константой или чем-то более разумным, учитывая ваши знания a и b.

Вы должны быть осторожны, потому что это не позволит идентифицировать a + b. Мы используем эту конструкцию как иерархическую априорную для симплексов:

parameters {
  real<lower = 1> kappa;
  real<lower = 0, upper = 1> phi;
  vector<lower = 0, upper = 1>[K] theta;

model {
  kappa ~ pareto(1, 1.5);  // power law prior
  phi ~ beta(a, b);  // choose your prior for theta
  theta ~ beta(kappa * phi, kappa * (1 - phi));  // vectorized

В моем тематическом исследовании Стэна повторных бинарных испытаний есть расширенный пример, доступный на странице тематических исследований на веб-сайте Стэна (каталог тематических исследований в настоящее время находится по ссылке на документацию на вкладке пользователей).

person Bob Carpenter    schedule 15.09.2017

arrow_upward
0
arrow_downward

Следуя предложениям в комментариях, я не уверен, что буду следовать этому подходу, но для справки я подумал, что по крайней мере опубликую ответ на свой вопрос о том, как это можно сделать в Стэне.

После некоторого расспроса о Stan Discourses и дальнейшего исследования я обнаружил, что решение состояло в том, чтобы установить собственное распределение плотности и использовать синтаксис target +=. Таким образом, эквивалент для Стэна примера для pymc будет:

parameters {
  real<lower=0> a;
  real<lower=0> b;
  real<lower=0,upper=1> p;
  ...
}

model {
  target += log((a + b)^-2.5);

  p ~ beta(a,b)
  ...
}
person Ian_Fin    schedule 13.09.2017