Как я могу (численно) приблизительные значения для бета-распределения с большой альфа и бета

Существует ли численно устойчивый способ расчета значений бета-распределения для большого целого числа альфа, бета (например, альфа, бета> 1000000)?

На самом деле, мне нужен только 99% доверительный интервал для режима, если это как-то облегчает проблему.

Добавить : Извините, мой вопрос был не так четко сформулирован, как я думал. Я хочу сделать следующее: у меня есть машина, которая проверяет продукты на ленточном конвейере. Некоторая часть этих продуктов отбраковывается машиной. Теперь, если оператор машины меняет какую-либо настройку проверки, я хочу показать ему / ей предполагаемый уровень брака и некоторые подсказки о том, насколько надежна текущая оценка.

Поэтому я подумал, что я рассматриваю фактическую частоту отклонений как случайную величину X и вычисляю распределение вероятностей для этой случайной величины на основе количества отклоненных объектов N и принятых объектов M. Если я предполагаю равномерное предварительное распределение для X, это бета-распределение в зависимости от N и M. Я могу либо отобразить это распределение непосредственно для пользователя, либо найти интервал [l, r], чтобы фактическая частота брака находилась в этом интервале с p> = 0,99 (используя терминологию шаббычефа), и отобразить это интервал. Для малых M, N (т. Е. Сразу после изменения параметра) я могу рассчитать распределение напрямую и приблизить интервал [l, r]. Но для больших M, N этот наивный подход приводит к ошибкам недостаточного значения, поскольку x ^ N * (1-x) ^ M слишком мало, чтобы его можно было представить как число с плавающей запятой двойной точности.

Я полагаю, что лучше всего использовать мое наивное бета-распределение для малых M, N и перейти к нормальному распределению с тем же средним и дисперсией, как только M, N превысит некоторый порог. Имеет ли это смысл?

Нормальное приближение работает очень хорошо, особенно в хвостах. Используйте среднее значение и дисперсию . Например, абсолютная относительная ошибка в вероятности хвоста в сложной ситуации (где может иметь место асимметрия), такой как достигает пика около и составляет менее когда вы более 1 SD от среднего. (Это не потому, что бета очень велика: при абсолютные относительные ошибки ограниченыα $\alpha/(\alpha+\beta)$ α=106,β=1080,000260,00006α=β=106$\frac{\alpha\beta}{(\alpha+\beta)^{2} (1+\alpha+\beta)}$$\alpha = 10^6, \beta = 10^8$$0.00026$$0.00006$$\alpha = \beta = 10^6$$0.0000001$.) Таким образом, это приближение отлично подходит практически для любых целей, включающих интервалы 99%.

В свете правок этого вопроса, обратите внимание, что бета-интегралы не вычисляются путем фактической интеграции подынтегральной функции: конечно, вы получите недочеты (хотя они на самом деле не имеют значения, поскольку они не вносят заметного вклада в интеграл) , Существует множество способов вычисления интеграла или его аппроксимации, как описано в Johnson & Kotz (Распределения в статистике). Онлайн-калькулятор можно найти по адресу http://www.danielsoper.com/statcalc/calc37.aspx . Вам действительно нужно обратное значение этого интеграла. Некоторые методы вычисления обратного описаны на сайте Mathematica по адресу http://functions.wolfram.com/GammaBetaErf/InverseBetaRegularized/., Код предоставлен в Числовых Рецептах (www.nr.com). Очень хороший онлайн калькулятор - это сайт Wolfram Alpha (www.wolframalpha.com): введите inverse beta regularized (.005, 1000000, 1000001)для левой конечной точки и inverse beta regularized (.995, 1000000, 1000001)для правой конечной точки ( , интервал 99%).$\alpha=1000000, \beta=1000001$

Быстрый графический эксперимент показывает, что бета-распределение очень похоже на нормальное распределение, когда альфа и бета очень велики. Погуглив «нормальный бета-предел распространения», я нашел http://nrich.maths.org/discus/messages/117730/143065.html?1200700623 , который дает «доказательство» ручной работы.

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

$[l,r]$$l$$r$$[l,r]$$\alpha, \beta$ $l$$r$ как отдельные цифры, так что этот маршрут может быть достаточно хорошим.

$p$$p$$log(p/(1-p))$$min(\alpha,\beta) > 100$

Например

f <- function(n, a, b) {
    p <- rbeta(n, a, b)
    lor <- log(p/(1-p))
    ks.test(lor, 'pnorm', mean(lor), sd(lor))$p.value
}
summary(replicate(50, f(10000, 100, 1000000)))

как правило, производит вывод, как

резюме (копия (50, f (10000, 100, 1000000))) Мин. 1 кв. Медиана Среднее 3 кв. Максимум. 0,01205 0,10870 0,18680 0,24810 0,36170 0,68730

т.е. типичные значения р составляют около 0,2.

$\alpha=100, \beta=100000$

$p$

f2 <- function(n, a, b) {
    p <- rbeta(n, a, b)
    ks.test(p, 'pnorm', mean(p), sd(p))$p.value
}
summary(replicate(50, f2(10000, 100, 1000000)))

производит что-то вроде

summary(replicate(50, f2(10000, 100, 1000000)))
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
2.462e-05 3.156e-03 7.614e-03 1.780e-02 1.699e-02 2.280e-01 

с типичными значениями р около 0,01

Функция R qqnormтакже дает полезную визуализацию, создавая очень прямолинейный график для распределения логарифмических шансов, указывающий приблизительную нормальность, распределение переменной бета dsitribute создает отличительную кривую, указывающую на ненормальность

$\alpha,\beta$