2

我正在使用 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)

这是我的功能的代码

def log_polya(Z,alpha):
    """
    Z is a vector of counts
    https://en.wikipedia.org/wiki/Dirichlet-multinomial_distribution
    http://mimno.infosci.cornell.edu/info6150/exercises/polya.pdf
    """
    if not isinstance(alpha,np.ndarray):
        alpha = np.array(alpha)
    if not isinstance(Z,np.ndarray):
        Z = np.array(Z)
    #Concentration Parameter
    A = sum(alpha)
    #Number of Datapoints
    N = sum(Z)
    return gammaln(A) - gammaln(N+A) + sum(gammaln(Z+alpha) - gammaln(alpha))

def llike(data,alpha):
    return sum(data.apply(log_polya,1,alpha=alpha))

def log_polya_derivative(Z,alpha):
    if not isinstance(alpha,np.ndarray):
        alpha = np.array(alpha)
    if not isinstance(Z,np.ndarray):
        Z = np.array(Z)
    if 0. in Z+alpha:
        Warning("invalid prior parameter,nans should be produced")
    #Concentration Parameter
    A = sum(alpha)
    #Number of Datapoints
    N = sum(Z)
    K = len(Z)
    return np.array([psi(A) - psi(N+A) + psi(Z[i]+alpha[i]) - psi(alpha[i]) for i in xrange(K)])


def gradient_llike(data,alpha):
     return np.array(data.apply(log_polya_derivative,1,alpha=alpha).sum(0))

更新：仍然对此感到好奇，但对于那些对此问题的工作实现感兴趣的人，以下用于实现 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
4

0 回答 0