Вычислите бета-функцию с отрицательными аргументами
Бета-функция (не бета-распределение) beta(a,b) может быть вычислена, когда a или b или оба отрицательны. Однако в реализации R разрешены только неотрицательные аргументы.
Есть ли способ обойти это?
Распространенное определение beta(a,b) — gamma(a)*gamma(b)/gamma(a+b); следовательно, мы могли бы использовать lgamma(), чтобы избежать переполнения:
exp(lgamma(a)+lgamma(b)-lgamma(a+b))
но для некоторых a или некоторых b, скажем, когда a <- -2.5; b <- 3, решение lgamma(a) должно быть комплексным числом, но R возвращает только действительную часть...
@CarlWitthoft: меня просят узнать все функции из неуказанных источников? Как я узнаю, что исчерпал все возможные библиотеки? Это невозможный подход.
На чем основано ваше утверждение, что lgamma(a) should be a complex number для a <- -2.5?
@user2554330, gsl::lngamma_complex(-2.5) возвращается -0.05624372-3.141593i
если бы я подумал об этом еще немного, я бы, вероятно, понял, почему мнимая часть lngamma(-2.5) равна -pi; меня также немного беспокоит то, что gsl::lngamma() с радостью возвращает результаты для отрицательных целочисленных значений x, где есть полюс...
gamma(-2.5) отрицательно, но документально подтверждено, что lgamma(-2.5) представляет собой логарифм абсолютного значения gamma(-2.5).
Ответы 1
pracma::gammaz или gsl::lngamma_complex, похоже, работают (и согласен с Wolfram Alpha).
## convert values to complex (not sure if this is necessary)
a <- -2.5+0i; b <- 3+0i
library(pracma)
Re(gammaz(a)*gammaz(b)/(gammaz(a+b))) ## -1.0667
Лучше (поскольку, как вы заметили, это позволяет избежать проблем с переполнением, когда a и b большие):
library(gsl)
exp((lngamma_complex(a)+lngamma_complex(b))-lngamma_complex(a+b))
## [1] -1.066667-1.30629e-16i
Для меня не очевидно, когда бета-функция имеет действительное значение; если вы были уверены, что результат будет реальным, вы можете использовать Re(), чтобы отбросить мнимую часть результата. В противном случае этот ответ может пригодиться для отбрасывания фазза с плавающей запятой в мнимую часть...
Работает как шарм. Спасибо.
Этот вопрос, вероятно, будет закрыт, если, по крайней мере, вы не перечислите каждый пакет R в CRAN (и биокондукторе), который предоставляет какую-то бета-функцию, и не покажете, что ни один из них не поддерживает ваши потребности.