n>1000的β二项似然计算

Question

在计算beta二项式可能性时，我正在努力解决一个数值精度问题。我的目标是估计某个数字d的概率y mod 10=d，假设y是二项式的（n，p），p是Beta（a，b）。我想给大n一个快速的解决方案，我的意思是至少1000个。有一件事似乎给了我合理的答案，那就是使用模拟。在

def npliketest_exact(n,digit,a,b):
    #draw 1000 values of p 
    probs = np.array(beta.rvs(a,b,size=1000))
    #create an array of numbers whose last digit is digit 
    digits = np.arange(digit,n+1,10)
    #create a function that calculates the pmf at x given p
    exact_func = lambda x,p: binom(n,p).pmf(x)
    #given p, the likeklihood of last digit "digit" is the sum over all entries in digits
    likelihood = lambda p: exact_func(digits,p).sum()
    #return the average of that likelihood over all the draws
    return np.vectorize(likelihood)(probs).mean()

np.random.seed(1)
print npliketest_exact(1000,9,1,1) #0.0992310195195

这也许没问题，但我担心这个策略的精确性。特别是如果有更好/更精确的方法来计算，我很想知道怎么做。在

我已经开始尝试用对数似然法来得出这个答案，但即使这样，我也遇到了数值稳定性问题。在

^{pr2}$

由于β1,1的平均值为0.5，从n=1000的β二项式中得到y=500的概率应该比得到9的概率大得多，但是上面的计算显示了一个可疑的常量值。在

我尝试过的另一件事，在stackoverflow的其他地方被建议来处理这个问题，就是使用一些聪明的技巧来支持数值稳定性，这些显然隐藏在scipy的betaln公式中。在

def binomln(n, k): #log of the binomial coefficient 
    # Assumes binom(n, k) >= 0
    return -betaln(1 + n - k, 1 + k) - log(n + 1)

def log_betabinom_exact(n,k,a,b):
    return binomln(n,k) + betaln(k+a,n-k+b) - betaln(a,b)

print exp(log_betabinom_exact(1000,9,1,1)) #.000999000999001
print exp(log_betabinom_exact(1000,500,1,1)) #0.000999000999001

同样的可疑常数。如有任何建议，将不胜感激。在这方面使用sympy会有什么帮助吗？在

****跟进

抱歉，伙计们，我犯了个愚蠢的错误，Beta（1，1）是统一的，所以我得到的结果是有意义的。尝试不同的参数会使不同的k值看起来不一样