Polya分布的极大似然估计

Question

我正在使用scipy为Polya发行版编写MLE。Nelder-Mead方法是有效的，但是当运行BFGS时，我得到了一个“由于精度损失而不一定达到期望的误差”的错误。Nelder-Mead方法似乎对我的需求来说太慢了（我有很多相当大的数据，比如1000个表，在某些情况下甚至可以达到10x10000）。我尝试过使用check_grad函数，结果在下面的例子中很小（顺序为10^-2），所以我不确定这是否意味着对数似然的梯度中存在缺陷，或者可能性只是非常强的峰值。不管怎么说，我一直在努力地看我的代码，我看不出问题所在。下面是一些重新创建问题的示例代码

#setup some data
from numpy.random import dirichlet, multinomial
from scipy.optimize import check_grad

alpha = [10,30,50]
p = pd.DataFrame(dirichlet(alpha,200))
data = p.apply(lambda x: multinomial(500,x),1)
a = np.array(data.mean(0))

#optimize
result = minimize(lambda a: -1*llike(data,exp(a)),
    x0=np.log(a),
    method='Nelder-Mead')
x0=result.x
result = minimize(lambda a: -1*llike(data,exp(a)),
    x0=x0,
    jac=lambda a: -1*gradient_llike(data,np.exp(a)),
    method='BFGS')
exp(result.x) #should be close to alpha
#uhoh, let's check that this is right.
check_grad(func=lambda a: -1*llike(data,a),grad=lambda a: -1*gradient_llike(data,a),x0=alpha)

这是我函数的代码

^{pr2}$

---

更新：仍然对此感到好奇，但是对于那些对这个问题的有效实现感兴趣的人来说，下面实现Minka不动点算法的代码似乎工作得很好（即，快速恢复接近真正dirichlet参数的值）。在

def minka_mle_polya(data):
    """
    http://research.microsoft.com/en-us/um/people/minka/papers/dirichlet/minka-dirichlet.pdf
    """
    data = np.array(data)
    K = np.shape(data)[1]
    alpha = np.array(data.mean(0))
    alpha_new = np.ndarray((K))
    precision = 10
    while precision > 10**-5:
        for k in range(K):
            A = sum(alpha)
            N = data.sum(1)
            numerator = sum(
                    psi(data[:,k]+alpha[k])-psi(alpha[k])
                    )
            denominator = sum(
                psi(N+A)-psi(A)
                )
            alpha_new[k] = alpha[k]*numerator/denominator
        precision = sum(abs(alpha_new - alpha))
        alpha_old = np.array(alpha)
        alpha = np.array(alpha_new)
        print "Gap", precision