Compute E_{q(z) q(x)} [log p(z) + log p(x | z) + log p(y | x, z)]
1234
def expected_log_joint_probability(self):
        
        # E_{q(z)}[log p(z)]
        from pyslds.util import expected_hmm_logprob
        elp = expected_hmm_logprob(
            self.pi_0, self.trans_matrix,
            (self.expected_states, self.expected_transcounts, self._normalizer))

        # E_{q(x)}[log p(y, x | z)]  is given by aBl
        # To get E_{q(x)}[ aBl ] we multiply and sum
        elp += np.sum(self.expected_states * self.vbem_aBl)
        return elp
339408
function(q,a,b,x,s,t){
        return __add(__rol(__add(__add(a,q),__add(x,t)),s),b);
    }
823265
def _get_cct(x, y, z):
        
        x_e = 0.3320
        y_e = 0.1858

        n = ((x / (x + z + z)) - x_e) / ((y / (x + z + z)) - y_e)

        a_0 = -949.86315
        a_1 = 6253.80338
        a_2 = 28.70599
        a_3 = 0.00004

        t_1 = 0.92159
        t_2 = 0.20039
        t_3 = 0.07125

        cct = a_0 + a_1 * numpy.exp(-n / t_1) + a_2 * numpy.exp(-n / t_2) + a_3 * numpy.exp(-n / t_3)
        return cct
138527
def log_prior(self):
        """"""
        for p, b in zip(self.parameter_vector, self.parameter_bounds):
            if b[0] is not None and p < b[0]:
                return -np.inf
            if b[1] is not None and p > b[1]:
                return -np.inf
        return 0.0
582471
def mp2q(p, q):
    
    p, q = flatten(p), flatten(q)
    entropy_dist = 1 / len(p)
    return sum(entropy_dist * np.nan_to_num((p ** 2) / q * np.log(p / q)))
262848
def _expectation(p, mean1, none1, mean2, none2, nghp=None):
    
    with params_as_tensors_for(mean1, mean2):
        e_xxt = p.cov + (p.mu[:, :, None] * p.mu[:, None, :])  # NxDxD
        e_A1t_xxt_A2 = tf.einsum("iq,nij,jz->nqz", mean1.A, e_xxt, mean2.A)  # NxQ1xQ2
        e_A1t_x_b2t = tf.einsum("iq,ni,z->nqz", mean1.A, p.mu, mean2.b)  # NxQ1xQ2
        e_b1_xt_A2 = tf.einsum("q,ni,iz->nqz", mean1.b, p.mu, mean2.A)  # NxQ1xQ2
        e_b1_b2t = mean1.b[:, None] * mean2.b[None, :]  # Q1xQ2

        return e_A1t_xxt_A2 + e_A1t_x_b2t + e_b1_xt_A2 + e_b1_b2t
1212160
private double logL(int k, long n, double x) {
        if (x == 0.0) x = 0.01;
        if (x == 1.0) x = 0.99;
        return k * Math.log(x) + (n-k) * Math.log(1-x);
    }
1212095
def expected_hmm_logprob(pi_0, trans_matrix, stats):
    
    E_z, sum_E_ztztp1T, _ = stats
    T, K = E_z.shape
    assert sum_E_ztztp1T.shape == (K, K)

    out = 0
    out += np.dot(E_z[0], np.log(pi_0))
    out += np.sum(sum_E_ztztp1T * np.log(trans_matrix))
    return out
794567
private static function logBeta($p, $q)
    {
        if ($p != self::$logBetaCacheP || $q != self::$logBetaCacheQ) {
            self::$logBetaCacheP = $p;
            self::$logBetaCacheQ = $q;
            if (($p <= 0.0) || ($q <= 0.0) || (($p + $q) > LOG_GAMMA_X_MAX_VALUE)) {
                self::$logBetaCacheResult = 0.0;
            } else {
                self::$logBetaCacheResult = self::logGamma($p) + self::logGamma($q) - self::logGamma($p + $q);
            }
        }
        return self::$logBetaCacheResult;
    }
1175213
def LogBinomialCoef(n, k):
    
    return n * log(n) - k * log(k) - (n - k) * log(n - k)
610643
def log_joint(self, x):
        
        T, D = self.T, self.D_latent
        assert x.shape == (T, D)

        # prior log p(x) -- quadratic terms
        J_diag, J_upper_diag = self.sparse_J_prior
        lp = -0.5 * np.sum(x * symm_block_tridiag_matmul(J_diag, J_upper_diag, x))

        # prior log p(x) -- linear terms
        _, h_init, log_Z_init = self.info_init_params
        _, _, _, h1, h2, log_Z_dyn = self.info_dynamics_params
        lp += x[0].dot(h_init)
        lp += np.sum(x[:-1] * h1)
        lp += np.sum(x[1:] * h2)

        # prior log p(x) -- normalization constants
        lp += log_Z_init
        lp += np.sum(log_Z_dyn)

        # likelihood log p(y | x)
        lp += self.log_conditional_likelihood(x)

        return lp
382572
function ecp_proj(R, N, p, l) {
	var x = N[0], y = N[1], z = N[2];
	var xr = R[0], yr = R[1], zr = R[2];
	var t0 = _ecp_t0, t1 = _ecp_t1, t2 = _ecp_t2;
	fp_mul(t2, t0, fp_sqr(t1, fp_inv(t0, z, p, l), p, l), p, l);
	fp_mul(xr, x, t1, p, l);
	fp_mul(yr, y, t2, p, l);
	mpn_zero(zr, l), zr[0] = 1;
	return R;
}
577702
def log(self):
        
        
        # Init
        norm = self.norm()
        vecNorm = self.x**2 + self.y**2 + self.z**2
        tmp = self.w / norm
        q = Quaternion()
        
        # Calculate
        q.w = np.log(norm)
        q.x = np.log(norm) * self.x * np.arccos(tmp) / vecNorm
        q.y = np.log(norm) * self.y * np.arccos(tmp) / vecNorm
        q.z = np.log(norm) * self.z * np.arccos(tmp) / vecNorm
        
        return q
494766
def _u0Eq(logu,delta,pot,E,Lz22):
    """"""
    u= numpy.exp(logu)
    sinh2u= numpy.sinh(u)**2.
    cosh2u= numpy.cosh(u)**2.
    dU= cosh2u*actionAngleStaeckel.potentialStaeckel(u,numpy.pi/2.,pot,delta)
    return -(E*sinh2u-dU-Lz22/delta**2./sinh2u)
693838
def qloguniform(low, high, q, random_state):
    '''
    
    '''
    return np.round(loguniform(low, high, random_state) / q) * q
892326
def grad_log_q(self,z):
                
        param_count = 0
        grad = np.zeros((np.sum(self.approx_param_no),self.sims))
        for core_param in range(len(self.q)):
            for approx_param in range(self.q[core_param].param_no):
                grad[param_count] = self.q[core_param].vi_score(z[core_param],approx_param)        
                param_count += 1
        return grad
352889
function bup(a, b, x, y, n) {
      var logSeries;

      logSeries = Math.log(series(function(i, v) {
            if (i === 0) { return 1; }
            return v * x * (1 + (b - 1) / (a + i));
         }, n));
      return a * Math.log(x) + b * Math.log(y) - lbeta(a, b) - Math.log(a) +
         logSeries;
   }