estimator.py 文件源码

def _calc_alpha_beta(log_d, alpha0, log_beta0, gamma, log_psi):
    log_beta = logsumexp(log_psi + log_d[:, :, :, np.newaxis], axis=0)
    log_beta = log_normalize(log_beta, axis=1)

    log_smooth = np.log(10)
    alpha = alpha0
    N = gamma.shape[0]
    zero = 1e-30

    gamma_digamma_sum = digamma(gamma.sum(axis=1))[:, np.newaxis]
    g_offset = (digamma(gamma) - gamma_digamma_sum).sum(axis=0) / N
    # using log
    def next_alpha(alpha):
        das = digamma(alpha.sum())
        g = alpha * N * (das - digamma(alpha) + g_offset)
        h = alpha * N * (das + g_offset)
        z = N * das
        x = (alpha * g / h).sum()
        w = (alpha ** 2 / h).sum()
        return np.exp(np.log(alpha) - (g - x * alpha / (1/z + w)) / h)

    return (alpha, log_beta)