代码之家 › 专栏 › 技术社区 › HJA24

用scipy.optimize和log似然寻找β二项式分布的α和β

binomial-theorem beta distribution scipy python

13

HJA24 · 技术社区 · 7 年前

分布是β二项式if 磷 ,成功的概率,在二项式分布中,具有β分布和形状参数 α±gt;0 和 α& gt;0 . 形状参数定义成功的概率。我想找出 γ± 和 γ 从β二项式分布的角度来最好地描述我的数据。我的数据集 players 包含有关点击次数的数据( H ,at bats的数目( 抗体 )以及转换( H/ab )很多棒球运动员。我在Juliend回答的帮助下估计了PDF。 Beta Binomial Function in 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()

现在,我想写一个最小化的函数 洛格利贝塔比诺姆

 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)

有人能告诉我我做错了什么吗?非常感谢您的帮助!!

1 回复 | 直到 7 年前

1

5

Davide Fiocco 7 年前

要注意的一件事是 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)

然后,您可以使用

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)

这将产生合理的结果。

你可以检查 How to properly fit a beta distribution in python? 如果你想进一步修改你的代码的话。