Как быстро сэмплировать X, если exp (X) ~ Gamma?

У меня есть простая проблема выборки, где мой внутренний цикл выглядит так:

v = sample_gamma(k, a)

где sample_gammaобразцы из гамма-распределения, чтобы сформировать образец Дирихле.

Это работает хорошо, но для некоторых значений k / a некоторые из последующих вычислений теряются.

Я адаптировал его для использования переменных пространства журнала:

v = log(sample_gamma(k, a))

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

Есть ли способ напрямую сэмплировать без использования медленных функций, таких как ? Я попробовал поискать в Google, но я даже не знаю, имеет ли этот дистрибутив общее имя (log-gamma?). $X, \exp(X) \sim \text{Gamma}$ $\log()$

sampling gamma-distribution luispedro
источник

Все, что вам нужно сделать, это разделить каждую гамму вариации на их сумму. Как же тогда происходит недосып? И как логарифм решает эту проблему (вы не можете вычислить сумму, не возводя обратно возведение в степень в любом случае)?

whuber

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

Луиспедро

Добавление журналов математически некорректно: оно соответствует умножению гамм, а не их добавлению. Да, вы можете получить рабочие результаты, но у них точно не будет дистрибутива Dirichlet! Опять же, какова природа исходного недостатка и какие расчеты вы делаете, когда это происходит? С какими фактическими ценностями вы работаете?

whuber

@whuber Я мог бы слишком упростить описание. Я делаю все, что я {т = гамма (а, б); сумма + = т; d [i] = log (t)}; logsum = log (сумма); forall i {d [i] - = logsum; }. Раньше это было недостаточно, если a было очень маленьким.

Луиспедро

Понял: для около 0 вы будете в беде, несмотря ни на что. Интересная проблема!

α

$\alpha$

whuber

Ответы:

Рассмотрим небольшой параметр формы около 0, такой как . В диапазоне между 0 и , приблизительно , так что гамма PDF приблизительно . Это может быть интегрировано в приблизительный CDF, . Обращая его, мы видим степень : огромный показатель. Для это вызывает некоторую вероятность недостаточного значения (значение двойной точности меньше , более или менее). Вот график вероятности недосыпания в зависимости от логарифма по основанию-десяти $\alpha$ $\alpha = 1/100$ $\alpha$ $e^{-\alpha}$ $1$ $x^{\alpha-1}dx / \Gamma(\alpha)$ $F_\alpha(x) = \frac{x^\alpha}{\alpha \Gamma(\alpha)}$ $1/\alpha$ $\alpha = 1/100$ $10^{-300}$ $\alpha$ :

введите описание изображения здесь

Одно из решений состоит в том, чтобы использовать это приближение для генерации логарифмических (гамма-вариаций): по сути, попробуйте сгенерировать гамма-вариацию и, если она слишком мала, сгенерировать ее логарифм из этого приблизительного распределения мощности (как показано ниже). (Повторяйте это до тех пор, пока журнал не окажется в пределах диапазона недостаточного заполнения, чтобы он являлся допустимой заменой исходной переменной пониженного значения.) Для вычисления Дирихле вычтите максимум всех логарифмов из каждого значения журнала: это неявно изменяет масштаб всех гамма изменяется, поэтому она не влияет на значения Дирихле. Обрабатывайте любой получающийся журнал, который слишком мал (скажем, меньше -100), как журнал с истинным нулем. Экспонировать другие журналы. Теперь вы можете продолжить без потери.

Это займет еще больше времени, чем раньше, но по крайней мере это сработает!

Чтобы сгенерировать приблизительный журнал Gamma варьировать с параметром формы , предварительно вычислите . Это легко, потому что существуют алгоритмы для непосредственного вычисления значений log Gamma . Создайте случайное число с плавающей точкой между 0 и 1, возьмите его логарифм, разделите на и добавьте к нему $\alpha$ $C = \log(\Gamma(\alpha)) + \log(\alpha)$ $\alpha$ $C$

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

редактировать

В другом ответе OP описывает метод, в котором степень равномерной переменной (переменная ) умножается на переменную . Это работает, потому что pdf совместного распределения этих двух переменных равен . Чтобы найти pdf из мы подставляем , делим на хакобейский и интегрируем из . Интеграл должен находиться в диапазоне от до потому что , откуда $1/\alpha$ $B(\alpha)$ $\Gamma(\alpha+1)$ $\left(\alpha x^{\alpha-1}\right) \left(y^{\alpha}e^{-y}dy/\Gamma(\alpha+1)\right)$ $z = xy$ $y \to z/x$ $x$ $x$ $z$ $\infty$ $0 \le y \le 1$

pdf (z) = \frac{α}{Γ (α + 1)} \int_{z}^{\infty} (x^{α} / x) e^{- x} (z / x)^{α - 1} d x d z = \frac{1}{Γ (α)} z^{α - 1} e^{- z} d z,

$\text{pdf}(z)=\frac{\alpha}{\Gamma(\alpha+1)}\int_z^{\infty}{\left(x^\alpha / x\right) e^{-x} (z/x)^{\alpha-1} dx} dz = \frac{1}{\Gamma(\alpha)}z^{\alpha-1}e^{-z}dz,$

который является pdf дистрибутива . $\Gamma(\alpha)$

Все дело в том, что когда , значение, полученное из , вряд ли будет недооценено, и, суммируя его log и раз логарифм независимой равномерной переменной, мы будет иметь журнал изменений . Журнал, скорее всего, будет очень отрицательным, но мы обойдем конструкцию его антилога, который будет потерян в представлении с плавающей запятой. $0 \lt \alpha \lt 1$ $\Gamma(\alpha+1)$ $1/\alpha$ $\Gamma(\alpha)$

Whuber
источник

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

\frac{Γ (α)}{Γ (α) + Γ (1)} \sim B e t a (α, 1)

$\frac{\Gamma(\alpha)}{\Gamma(\alpha)+\Gamma(1)}\sim Beta(\alpha,1)$

Γ (α) + Γ (1) \sim Γ (α + 1)

$\Gamma(\alpha)+\Gamma(1)\sim \Gamma(\alpha+1)$

α \approx 0

$\alpha\approx 0$

y \sim e x p o (1)

$y\sim expo(1)$

- \log (u)

$-\log(u)$

Γ (α + 1)

$\Gamma(\alpha+1)$

вероятностная

Я отвечаю на свой вопрос, но нашел довольно хорошее решение, даже если не до конца его понимаю. Глядя на код из Научной библиотеки GNU, мы рассмотрим примеры гамма-переменных ( rэто генератор случайных чисел, aэто и is ): $\alpha$ b $\beta$

  if (a < 1)
    {
      double u = gsl_rng_uniform_pos (r);
      return gsl_ran_gamma (r, 1.0 + a, b) * pow (u, 1.0 / a);
   }

gsl_ran_gammaявляется функцией, которая возвращает гамма-случайную выборку (поэтому вышеупомянутый является рекурсивным вызовом), в то время как gsl_rng_uniform_posвозвращает равномерно распределенное число в ( для строго положительного, поскольку гарантированно не вернет 0.0). $(0,1)$ _pos

Поэтому я могу взять журнал последнего выражения и использовать

return log(gsl_ran_gamma(r, 1.0 + a, b)) + log(u)/a;

Чтобы получить то, что я хотел. У меня сейчас два log()звонка (но на один меньше pow()), но результат, вероятно, лучше. До того, как указал Уубер, у меня было что-то повышенное до , потенциально огромное количество. Теперь в logspace я умножаю на . Таким образом, это меньше шансов на снижение. $1/a$ $1/a$

luispedro
источник

Не могли бы вы объяснить, что делают gsl_rng_uniform_pos и gsl_ran_gamma? Я предполагаю, что первое возвращает равномерное случайное значение между 0 и r, а второе связано со значением гаммы (1 + a, b) - может быть, это неполная гамма? В целом это выглядит очень близко к предложенному мною приближению (за исключением того, что при его рассмотрении очевидно, что я забыл указать деление на часть, что очень важно!)

α

$\alpha$

whuber

Я отредактировал свой ответ, чтобы включить больше деталей сейчас.

Луиспедро

Спасибо: а что такое "р"? (Обратите внимание, что рекурсия ограничена: будет сделано не более одного рекурсивного вызова, потому что a> 0 подразумевает 1.0 + a> 1.)

whuber

r - генератор случайных чисел (откуда вы получаете случайные числа).

Луиспедро

Ах, это умно: произведение переменной и независимой переменной оказывается переменной . Я отредактировал свой ответ, чтобы он указал на ваше решение и объяснил, почему оно работает.

Γ (α + 1)

$\Gamma(\alpha+1)$

B (α, 1)

$B(\alpha,1)$

Γ (α)

$\Gamma(\alpha)$

whuber