Нахождение альфа и бета бета-биномиального распределения с scipy.optimize и log правдоподобия

Question

Распределение является бета-биномиальным, если p , вероятность успеха, в биномиальном распределении имеет бета-распределение с параметрами формы α> 0 и β> 0.Параметры формы определяют вероятность успеха.Я хочу найти значения для α и β , которые лучше всего описывают мои данные с точки зрения бета-биномиального распределения.Мой набор данных players состоит из данных о количестве обращений ( H ), количестве бит-бит ( AB ) и конверсии ( H / AB * 1017).*) из многих бейсболистов.Я оцениваю PDF с помощью ответа JulienD в Бета биноминальная функция в Python

from scipy.special import beta
from scipy.misc import comb

pdf = comb(n, k) * beta(k + a, n - k + b) / beta(a, b)

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

def loglike_betabinom(params, *args):
   """
   Negative log likelihood function for betabinomial distribution
   :param params: list for parameters to be fitted.
   :param args:  2-element array containing the sample data.
   :return: negative log-likelihood to be minimized.
   """

   a, b = params[0], params[1]
   k = args[0] # the conversion rate
   n = args[1] # the number of at-bats (AE)

   pdf = comb(n, k) * beta(k + a, n - k + b) / beta(a, b)

   return -1 * np.log(pdf).sum()   

Теперь я хочу написать функцию, которая минимизирует loglike_betabinom

 from scipy.optimize import minimize
 init_params = [1, 10]
 res = minimize(loglike_betabinom, x0=init_params,
                args=(players['H'] / players['AB'], players['AB']),
                bounds=bounds,
                method='L-BFGS-B',
                options={'disp': True, 'maxiter': 250})
 print(res.x)

Результат - [-6.04544138 2.03984464], что подразумевает, что α отрицателен, что невозможно.Я основал свой сценарий на следующем R-фрагменте.Они получают [101.359, 287.318] ..

 ll <- function(alpha, beta) { 
    x <- career_filtered$H
    total <- career_filtered$AB
    -sum(VGAM::dbetabinom.ab(x, total, alpha, beta, log=True))
 }

 m <- mle(ll, start = list(alpha = 1, beta = 10), 
 method = "L-BFGS-B", lower = c(0.0001, 0.1))

 ab <- coef(m)

Может кто-нибудь сказать мне, что я делаю неправильно?Помощь очень ценится!

Answer 1

Стоит обратить внимание на то, что comb(n, k) в вашем логарифмическом правдоподобии может не иметь численного поведения для значений n и k в вашем наборе данных.Вы можете проверить это, применив comb к вашим данным и посмотрите, появятся ли inf s.

Один из способов исправить ситуацию может состоять в том, чтобы переписать отрицательное логарифмическое правдоподобие, как предложено в https://stackoverflow.com/a/32355701/4240413,, то есть как функцию логарифмов гамма-функций, как в

from scipy.special import gammaln
import numpy as np

def loglike_betabinom(params, *args):

    a, b = params[0], params[1]
    k = args[0] # the OVERALL conversions
    n = args[1] # the number of at-bats (AE)

    logpdf = gammaln(n+1) + gammaln(k+a) + gammaln(n-k+b) + gammaln(a+b) - \
     (gammaln(k+1) + gammaln(n-k+1) + gammaln(a) + gammaln(b) + gammaln(n+a+b))

    return -np.sum(logpdf) 

Youзатем можно минимизировать логарифмическую вероятность с помощью

from scipy.optimize import minimize

init_params = [1, 10]
# note that I am putting 'H' in the args
res = minimize(loglike_betabinom, x0=init_params,
            args=(players['H'], players['AB']),
            method='L-BFGS-B', options={'disp': True, 'maxiter': 250})
print(res)

, и это должно дать разумные результаты.

Вы можете проверить Как правильно разместить бета-дистрибутив в python? длявдохновение, если вы хотите доработать свой код.