def _distn_form_heldout_log_likelihood(self, X, M=10):
        # assert len(self.data_list) == 1, "TODO: Support more than 1 data set"

        T, K = X.shape
        assert K == self.K
        kappa = kappa_vec(X)
        N = N_vec(X)

        # Compute the data-specific normalization constant from the
        # augmented multinomial distribution
        Z_mul = (gammaln(N + 1) - gammaln(X[:,:-1]+1) - gammaln(N-X[:,:-1]+1)).sum()
        Z_mul += (-N * np.log(2.)).sum()

        # Monte carlo integrate wrt omega ~ PG(N, 0)
        import pypolyagamma as ppg
        hlls = np.zeros(M)
        for m in xrange(M):
            # Sample omega using the emission distributions samplers
            omega = np.zeros(N.size)
            ppg.pgdrawvpar(self.emission_distn.ppgs,
                           N.ravel(), np.zeros(N.size),
                           omega)
            omega = omega.reshape((T, K-1))
            # valid = omega > 0
            omega = np.clip(omega, 1e-8, np.inf)

            # Compute the normalization constant for this omega
            Z_omg = 0.5 * (kappa**2/omega).sum()
            Z_omg += T * (K-1)/2. * np.log(2*np.pi)
            Z_omg += -0.5 * np.log(omega).sum()  # 1/2 log det of Omega_t^{-1}

            # clip omega = zero for message passing
            # omega[~valid] = 1e-32

            # Exactly integrate out the latent states z using message passing
            # The "data" is the normal potential from the
            states = MultinomialLDSStates(model=self, data=X)
            conditional_mean = kappa / omega - self.emission_distn.mu[None, :]
            # conditional_mean[~np.isfinite(conditional_mean)] = 0
            conditional_cov = np.zeros((T, K-1, K-1))
            for t in xrange(T):
                conditional_cov[t,:,:] = np.diag(1./omega[t,:])
            Z_lds = states.log_likelihood(conditional_mean, conditional_cov)

            # Sum them up to get the heldout log likelihood for this omega
            hlls[m] = Z_mul + Z_omg + Z_lds

        # Now take the log of the average to get the log likelihood
        hll = logsumexp(hlls) - np.log(M)

        # Use bootstrap to compute error bars
        samples = np.random.choice(hlls, size=(100, M), replace=True)
        hll_samples = logsumexp(samples, axis=1) - np.log(M)
        std_hll = hll_samples.std()

        return hll, std_hll
File to file
dataset/ETHPy150Open HIPS/pgmult/pgmult/lds.py/_MultinomialLDSBase._distn_form_heldout_log_likelihood--crossover with--dataset/ETHPy150Open mapbox/rasterio/tests/test_warp.py/test_warp_from_to_file_multi