的分配是β-二项式如果p,成功的概率,在二项式分布具有带形状参数α> 0和β> 0形状参数定义成功的概率β分布。我想找到的α和β从β-二项式分布的角度来看最能描述我的数据值。我的数据集players
大约命中(H)的数量,在蝙蝠(AB)和大量的棒球选手的转换(H / AB)的数量由数据。我估计PDF与JulienD的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()
现在,我想写最小化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)
谁能告诉我什么,我做错了什么?帮助非常感谢!
有一点要注意的是在您的登录可能性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?灵感,如果你想进一步返工你的代码。