vb.py 文件源码

def update_tau(self):
        sum_nu = numpy.sum(self._nu, axis=0)
        assert(len(sum_nu) == self._K)
        if self._finite_mode:
            assert(len(sum_nu) == self._K)
            self._tau[0, :] = self._alpha / self._K + sum_nu
            self._tau[1, :] = 1. + self._N - sum_nu
        else:
            N_minus_sum_nu = self._N - sum_nu
            for k in xrange(self._K):
                psi_tau = scipy.special.psi(self._tau)
                assert(psi_tau.shape == ( self._K,2))
                #assert(psi_tau.shape == (2, self._K))

                psi_sum_tau = scipy.special.psi(numpy.sum(self._tau, axis=0))
                assert(len(psi_sum_tau) == self._K)

                psi_tau0_cumsum = numpy.hstack([0, numpy.cumsum(psi_tau[0, :-1])])
                assert(len(psi_tau0_cumsum) == self._K)

                psi_sum_cumsum = numpy.cumsum(psi_sum_tau)
                assert(len(psi_sum_cumsum) == self._K)

                exponent = psi_tau[1, :] + psi_tau0_cumsum - psi_sum_cumsum
                unnormalized = numpy.exp(exponent - numpy.max(exponent))
                assert(len(unnormalized) == self._K)

                qs = numpy.zeros((self._K, self._K))
                for m in xrange(k, self._K):
                    qs[m, 0:m + 1] = unnormalized[0:m + 1] / numpy.sum(unnormalized[0:m + 1])

                self._tau[0, k] = numpy.sum(sum_nu[k:self._K]) + numpy.dot(N_minus_sum_nu[k + 1:self._K], numpy.sum(qs[k + 1:self._K, k + 1:self._K], axis=1)) + self._alpha
                self._tau[1, k] = numpy.dot(N_minus_sum_nu[k:self._K], qs[k:self._K, k]) + 1
        return