Posterior covariance
        $$
        K_{xx} - K_{xx}W_{xx}^{-1}K_{xx}
        W_{xx} := \texttt{Woodbury inv}
        $$
421304
def covariance(self):
        
        if self._covariance is None:
            # LiK, _ = dtrtrs(self.woodbury_chol, self._K, lower=1)
            self._covariance = (
            np.atleast_3d(self._K) - np.tensordot(np.dot(np.atleast_3d(self.woodbury_inv).T, self._K), self._K,
                                                  [1, 0]).T).squeeze()
            # self._covariance = self._K - self._K.dot(self.woodbury_inv).dot(self._K)
        return self._covariance
515222
def _forward_step(self, Y, T_X, Var_X, Var_dX, rho_e, sigma2_e, weight):
        
        # We currently only implement diagonal form
        # of covariance matrix for Var_X, Var_dX and T_X, which means
        # each dimension of X is independent and their innovation noise
        # are also independent. Note that log_p_data takes this assumption.
        if Var_X.ndim == 1:
            inv_Var_X = np.diag(1 / Var_X)
            half_log_det_Var_X = np.sum(np.log(Var_X)) / 2.0
            Var_X = np.diag(Var_X)
            # the marginal variance of X
        else:
            half_log_det_Var_X = self._half_log_det(Var_X)
            inv_Var_X = np.linalg.inv(Var_X)
        if Var_dX.ndim == 1:
            inv_Var_dX = np.diag(1 / Var_dX)
            half_log_det_Var_dX = np.sum(np.log(Var_dX)) / 2.0
            Var_dX = np.diag(Var_dX)
            # the marginal variance of Delta X (the change of X from
            # previous time point)
        else:
            inv_Var_dX = np.linalg.inv(Var_dX)
            half_log_det_Var_dX = self._half_log_det(Var_dX)
        if T_X.ndim == 1:
            T_X = np.diag(T_X)
            # Transfer function of X: the expected mean of X at t+1
            # time point is T_x * X
        [n_T, n_V] = np.shape(Y)
        # numbers of time points and voxels
        mu = [None] * n_T
        # posterior mean of X, conditioned on all data up till the current
        # time point
        Gamma_inv = [None] * n_T
        # inverse of poterior Gamma.
        mu_Gamma_inv = [None] * n_T
        # mu * inv(Gamma)
        log_p_data = - np.log(np.pi * 2) * (n_T * n_V) / 2 \
            - half_log_det_Var_X - np.sum(np.log(sigma2_e)) * n_T / 2.0\
            + np.sum(np.log(1 - rho_e**2)) / 2.0 - half_log_det_Var_dX \
            * (n_T - 1)
        # This is the term to be incremented by c_n at each time step.
        # We first add all the fixed terms to it.

        # The following are a few fixed terms.
        Lambda_0 = np.dot(T_X, np.dot(inv_Var_dX, T_X.T)) \
            + np.dot(weight * rho_e**2 / sigma2_e, weight.T)
        H = np.dot(inv_Var_dX, T_X.T) + np.dot(weight * rho_e / sigma2_e,
                                               weight.T)
        Lambda_1 = inv_Var_dX + np.dot(weight / sigma2_e, weight.T)

        Gamma_inv[0] = inv_Var_X + np.dot(
            weight * (1 - rho_e**2) / sigma2_e, weight.T)
        # We might not need this and only use linalg.solve for related terms.
        mu_Gamma_inv[0] = np.dot(
            Y[0, :] * (1 - rho_e**2) / sigma2_e, weight.T)
        mu[0] = np.linalg.solve(Gamma_inv[0], mu_Gamma_inv[0])
        log_p_data -= 0.5 * np.sum(Y[0, :]**2 * (1 - rho_e**2) / sigma2_e)
        # This is the term added for the first time point.

        deltaY = Y[1:, :] - rho_e * Y[:-1, :]
        deltaY_sigma2inv_rho_weightT = np.dot(
            deltaY / sigma2_e * rho_e, weight.T)
        for t in np.arange(1, n_T):
            Gamma_tilde_inv = Lambda_0 + Gamma_inv[t - 1]
            tmp = np.linalg.solve(Gamma_tilde_inv, H.T)
            Gamma_inv[t] = Lambda_1 - np.dot(H, tmp)
            mu_Gamma_inv[t] = np.dot(deltaY[t - 1, :] / sigma2_e, weight.T) \
                + np.dot(mu_Gamma_inv[t - 1]
                         - deltaY_sigma2inv_rho_weightT[t - 1, :], tmp)
            mu[t] = np.linalg.solve(Gamma_inv[t], mu_Gamma_inv[t])
            tmp2 = mu_Gamma_inv[t - 1] - deltaY_sigma2inv_rho_weightT[t - 1, :]
            log_p_data += -self._half_log_det(Gamma_tilde_inv) \
                + np.dot(tmp2, np.linalg.solve(Gamma_tilde_inv, tmp2)) / 2.0
        log_p_data += -self._half_log_det(Gamma_inv[-1]) \
            + np.dot(mu_Gamma_inv[-1], mu[-1]) / 2.0 \
            - np.sum(deltaY**2 / sigma2_e) / 2.0
        return mu, mu_Gamma_inv, Gamma_inv, log_p_data, Lambda_0, \
            Lambda_1, H, deltaY, deltaY_sigma2inv_rho_weightT
955951
def map_parameter_estimate(self):
        
        res = self.pst.res
        assert res is not None
        # build the prior expectation parameter vector
        prior_expt = self.pst.parameter_data.loc[:,["parval1"]].copy()
        islog = self.pst.parameter_data.partrans == "log"
        prior_expt.loc[islog] = prior_expt.loc[islog].apply(np.log10)
        prior_expt = Matrix.from_dataframe(prior_expt)
        prior_expt.col_names = ["prior_expt"]
        # build the residual vector
        res_vec = Matrix.from_dataframe(res.loc[:,["residual"]])

        # form the terms of Schur's complement
        b = self.parcov * self.jco.T
        c = ((self.jco * self.parcov * self.jco.T) + self.obscov).inv
        bc = Matrix((b * c).x, row_names=b.row_names, col_names=c.col_names)

        # calc posterior expectation
        upgrade = bc * res_vec
        upgrade.col_names = ["prior_expt"]
        post_expt = prior_expt + upgrade

        # post processing - back log transform
        post_expt = pd.DataFrame(data=post_expt.x,index=post_expt.row_names,
                                 columns=["post_expt"])
        post_expt.loc[:,"prior_expt"] = prior_expt.x.flatten()
        post_expt.loc[islog,:] = 10.0**post_expt.loc[islog,:]
        # unc_sum = self.get_parameter_summary()
        # post_expt.loc[:,"standard_deviation"] = unc_sum.post_var.apply(np.sqrt)
        post_expt.sort_index(inplace=True)
        return post_expt
394483
def stochastic_kl_builder(posterior, prior, sample):
  """"""
  return tf.subtract(
      tf.reduce_sum(posterior.log_prob(sample)),
      tf.reduce_sum(prior.log_prob(sample)))
393534
def covariance2(xs, ys)
      raise "Length xs and ys must be equal" unless xs.length == ys.length

      ev_x, ev_y = mean(xs), mean(ys)
      xs.zip(ys)
        .map{|x,y| (x-ev_x) * (y-ev_y)}
        .inject(0) {|sum, x| sum += x} / xs.length
    end
538517
def rwishart_cov(n, C):
    
    # return rwishart(n, np.linalg.inv(C))

    p = np.shape(C)[0]
    # Need cholesky decomposition of precision matrix C^-1?
    sig = np.linalg.cholesky(C)

    if n <= (p-1):
        raise ValueError('Wishart parameter n must be greater '
                         'than size of matrix.')

    norms = np.random.normal(size=(p * (p - 1)) // 2)
    chi_sqs = np.sqrt(np.random.chisquare(df=np.arange(n, n - p, -1)))
    A = flib.expand_triangular(chi_sqs, norms)

    flib.dtrmm_wrap(sig, A, side='L', uplo='L', transa='N', alpha=1.)
    w = np.asmatrix(np.dot(A, A.T))
    flib.symmetrize(w)
    return w
179612
def _mse_converged(self):
        

        mse = mean_squared_error(self.local_prior, self.local_posterior_,
                                 multioutput='uniform_average')
        if mse > self.threshold:
            return False, mse
        else:
            return True, mse
941429
def q_prior(q, m=1, gamma=0.3, qmin=0.1):
    
    if q < qmin or q > 1:
        return 0
    C = 1/(1/(gamma+1)*(1 - qmin**(gamma+1)))
    return C*q**gamma
1121405
private double joint(int k, int y[]) {
        double p = 1.0; //priors[k];
        for (int j = 0; j < y.length; j++) {
            if (j != k && y[j] == 1) {
                p *= Conditional[k][j];
            }
        }
        return p;
    }
125823
def MultivariateGaussian(mu: vertex_constructor_param_types, covariance: vertex_constructor_param_types, label: Optional[str]=None) -> Vertex:
    
    return Double(context.jvm_view().MultivariateGaussianVertex, label, cast_to_double_vertex(mu), cast_to_double_vertex(covariance))
131032
def _u_distance_covariance_sqr_naive(x, y, exponent=1):
    
    a = _u_distance_matrix(x, exponent=exponent)
    b = _u_distance_matrix(y, exponent=exponent)

    return u_product(a, b)
630441
def generate_covariance(ts, sigma, tau):

    
    """

    ndim = ts.shape[0]

    tis = ts[:, np.newaxis]
    tjs = ts[np.newaxis, :]

    return sigma*sigma*np.exp(-np.square(tis-tjs)/(2.0*tau*tau))
1190905
def lnprior(pars):
    
    # Limit norm and B to be positive
    logprob = (
        naima.uniform_prior(pars[0], 0.0, np.inf)
        + naima.uniform_prior(pars[1], -1, 5)
        + naima.uniform_prior(pars[3], 0, np.inf)
    )

    return logprob
1289088
def deconv2d(
    input_, output_shape, k_h, k_w, d_h, d_w, stddev=0.02, name="deconv2d"):
  """"""
  with tf.variable_scope(name):
    w = tf.get_variable(
        "w", [k_h, k_w, output_shape[-1], input_.get_shape()[-1]],
        initializer=tf.random_normal_initializer(stddev=stddev))
    deconv = tf.nn.conv2d_transpose(
        input_, w, output_shape=output_shape, strides=[1, d_h, d_w, 1])
    biases = tf.get_variable(
        "biases", [output_shape[-1]], initializer=tf.constant_initializer(0.0))
    return tf.reshape(tf.nn.bias_add(deconv, biases), deconv.get_shape())
611029
def decompose_covariance(c):
    

    # make it a kickass copy of the original
    c = _n.array(c)

    # first get the error vector
    e = []
    for n in range(0, len(c[0])): e.append(_n.sqrt(c[n][n]))

    # now cycle through the matrix, dividing by e[1]*e[2]
    for n in range(0, len(c[0])):
        for m in range(0, len(c[0])):
            c[n][m] = c[n][m] / (e[n]*e[m])

    return [_n.array(e), _n.array(c)]
1082753
def calc_K_xx_wo_variance(self,X):
        """"""
        K_xx = np.ones([X.shape[0],X.shape[0]]) #ones now as a product occurs over each dimension
        for i,x in enumerate(X):
            for j,x2 in enumerate(X):
                for il,l in enumerate(self.lengthscale):
                    idx = il*2 #each pair of input dimensions describe the limits on one actual dimension in the data
                    K_xx[i,j] *= self.k_xx(x[idx],x2[idx],x[idx+1],x2[idx+1],l)
        return K_xx
1319090
def noise_covariance(fit, dof=2, **kw):
    
    ev = fit.eigenvalues

    measurement_noise = ev[-1]/(fit.n-dof)
    return 4*ev*measurement_noise