This takes an error vector and a correlation matrix and assembles the covariance
10490
def assemble_covariance(error, correlation):
    

    covariance = []
    for n in range(0, len(error)):
        covariance.append([])
        for m in range(0, len(error)):
            covariance[n].append(correlation[n][m]*error[n]*error[m])
    return _n.array(covariance)
327156
def _refit(self):
        
        # temporary variables for manipulation
        mu_diff = np.empty_like(self.f.components[0].mu)
        sigma   = np.empty_like(self.f.components[0].sigma)
        mean    = np.empty_like(mu_diff)
        cov     = np.empty_like(sigma)

        for j, c in enumerate(self.g.components):
            # stop if inv_map is empty for j-th comp.
            if not self.inv_map[j]:
                self.g.weights[j] = 0.
                continue

            # (re-)initialize new mean/cov to zero
            mean[:] = 0.0
            cov[:] = 0.0

            # compute total weight and mean
            self.g.weights[j] = self.f.weights[self.inv_map[j]].sum()
            for i in self.inv_map[j]:
                mean += self.f.weights[i] * self.f.components[i].mu

            # rescale by total weight
            mean /= self.g.weights[j]

            # update covariance
            for i in self.inv_map[j]:
                # mu_diff = mu'_j - mu_i
                mu_diff[:] = mean
                mu_diff -= self.f.components[i].mu

                # sigma = (mu'_j - mu_i) (mu'_j - mu_i)^T
                sigma[:] = np.outer(mu_diff, mu_diff)

                # sigma += sigma_i
                sigma += self.f.components[i].sigma

                # multiply with alpha_i
                sigma *= self.f.weights[i]

                # sigma_j += alpha_i * (sigma_i + (mu'_j - mu_i) (mu'_j - mu_i)^T
                cov += sigma

            # 1 / beta_j
            cov /= self.g.weights[j]

            # update the Mixture
            c.update(mean, cov)
977202
def iterative_plane_errors(axes,covariance_matrix, **kwargs):
    
    sheet = kwargs.pop('sheet','upper')
    level = kwargs.pop('level',1)
    n = kwargs.pop('n',100)

    cov = N.sqrt(N.diagonal(covariance_matrix))
    u = N.linspace(0, 2*N.pi, n)

    scales = dict(upper=1,lower=-1,nominal=0)
    c1 = scales[sheet]
    c1 *= -1 # We assume upper hemisphere
    if axes[2,2] < 0:
        c1 *= -1

    def sdot(a,b):
        return sum([i*j for i,j in zip(a,b)])

    def step_func(a):
        e = [
            N.cos(a)*cov[0],
            N.sin(a)*cov[1],
            c1*cov[2]]
        d = [sdot(e,i)
            for i in axes.T]
        x,y,z = d[2],d[0],d[1]
        r = N.sqrt(x**2 + y**2 + z**2)
        lat = N.arcsin(z/r)
        lon = N.arctan2(y, x)
        return lon,lat

    # Get a bundle of vectors defining
    # a full rotation around the unit circle
    return N.array([step_func(i)
        for i in u])
582889
def compute_symm_block_tridiag_covariances(H_diag, H_upper_diag):
    
    T, D, _ = H_diag.shape
    assert H_diag.ndim == 3 and H_diag.shape[2] == D
    assert H_upper_diag.shape == (T - 1, D, D)

    J_init = J_11 = J_22 = np.zeros((D, D))
    h_init = h_1 = h_2 = np.zeros((D,))

    J_21 = np.swapaxes(H_upper_diag, -1, -2)
    J_node = H_diag
    h_node = np.zeros((T, D))

    _, _, sigmas, E_xt_xtp1 = \
        info_E_step(J_init, h_init, 0,
                    J_11, J_21, J_22, h_1, h_2, np.zeros((T-1)),
                    J_node, h_node, np.zeros(T))
    return sigmas, E_xt_xtp1
941332
def readCovarianceMatrixFile(cfile,readCov=True,readEig=True):
    "
    covFile = cfile+'.cov'
    evalFile = cfile+'.cov.eval'
    evecFile = cfile+'.cov.evec'

    RV = {}
    if readCov:
        assert os.path.exists(covFile), '%s is missing.'%covFile
        RV['K'] = SP.loadtxt(covFile)
    if readEig:
        assert os.path.exists(evalFile), '%s is missing.'%evalFile
        assert os.path.exists(evecFile), '%s is missing.'%evecFile
        RV['eval'] = SP.loadtxt(evalFile)
        RV['evec'] = SP.loadtxt(evecFile)

    return RV
1359199
public static Matrix calculateCovarianceMatrixOfCenteredData(Matrix data){
		Assume.assume(data.getRowDimension()>1, "Number of data samples is "+
				data.getRowDimension()+", but it has to be >1 to compute "+
				"covariances");
		int dimsNo = data.getColumnDimension();
		int samplesNo = data.getRowDimension();
		Matrix m = new Matrix(dimsNo, dimsNo);
		for(int r = 0; r < dimsNo; r++)
			for(int c = r; c < dimsNo; c++){
				double sum = 0;
				for(int i = 0; i < samplesNo; i++)
					 sum += data.get(i, r)*data.get(i, c);
				m.set(r, c, sum/(samplesNo-1));
			}
		for(int r = 0; r < dimsNo; r++)
			for(int c = 0; c < r; c++) m.set(r, c, m.get(c, r));
		return m;		
	}
}
1322561
def setCovariance(self,cov):
        
        chol = LA.cholesky(cov,lower=True)
        params = chol[sp.tril_indices(self.dim)]
        self.setParams(params)
770538
def _calc_covar_matrix(self, profile):
        
        corr = self._calc_corr(profile)
        std = self._calc_ln_std(profile)
        # Modify the standard deviation by the truncated norm scale
        std *= randnorm.scale

        var = std ** 2
        covar = corr * std[:-1] * std[1:]

        # Main diagonal is the variance
        mat = diags([covar, var, covar], [-1, 0, 1]).toarray()

        return mat
1361461
def eighCovarianceMatrix(cfile):
    
    # precompute eigenvalue decomposition
    K = np.loadtxt(cfile+'.cov')
    K+= 1e-4*sp.eye(K.shape[0])
    S,U = la.eigh(K); S=S[::-1]; U=U[:,::-1]
    np.savetxt(cfile+'.cov.eval',S,fmt='%.6f')
    np.savetxt(cfile+'.cov.evec',U,fmt='%.6f')
1252443
def perform_pca(A):
    
    # First subtract the mean
    M = (A-numpy.mean(A.T, axis=1)).T
    # Get eigenvectors and values of covariance matrix
    return numpy.linalg.eig(numpy.cov(M))
139995
public static void recomputeCovarianceMatrices(Relation<? extends NumberVector> relation, WritableDataStore<double[]> probClusterIGivenX, List<? extends EMClusterModel<?>> models, double prior) {
    final int k = models.size();
    boolean needsTwoPass = false;
    for(EMClusterModel<?> m : models) {
      m.beginEStep();
      needsTwoPass |= m.needsTwoPass();
    }
    // First pass, only for two-pass models.
    if(needsTwoPass) {
      for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
        double[] clusterProbabilities = probClusterIGivenX.get(iditer);
        NumberVector instance = relation.get(iditer);
        for(int i = 0; i < clusterProbabilities.length; i++) {
          final double prob = clusterProbabilities[i];
          if(prob > 1e-10) {
            models.get(i).firstPassE(instance, prob);
          }
        }
      }
      for(EMClusterModel<?> m : models) {
        m.finalizeFirstPassE();
      }
    }
    double[] wsum = new double[k];
    for(DBIDIter iditer = relation.iterDBIDs(); iditer.valid(); iditer.advance()) {
      double[] clusterProbabilities = probClusterIGivenX.get(iditer);
      NumberVector instance = relation.get(iditer);
      for(int i = 0; i < clusterProbabilities.length; i++) {
        final double prob = clusterProbabilities[i];
        if(prob > 1e-10) {
          models.get(i).updateE(instance, prob);
        }
        wsum[i] += prob;
      }
    }
    for(int i = 0; i < models.size(); i++) {
      // MLE / MAP
      final double weight = prior <= 0. ? wsum[i] / relation.size() : (wsum[i] + prior - 1) / (relation.size() + prior * k - k);
      models.get(i).finalizeEStep(weight, prior);
    }
  }
734187
private function computeCovariance()
    {
        $factorCallable = $this->coVarianceFactor;
        $this->factor = $factorCallable();

        // Cache covariance and inverse covariance of the data
        if (null === $this->_dataInvCov) {
            // original used the numpy.cov function.
            $this->_dataCovariance = pow(Statistics::stdev($this->dataset, true), 2);

            //$this->_dataInvCov = 1/ linalg.inv($this->_dataCovariance)
            $this->_dataInvCov = 1 / $this->_dataCovariance;
        }

        $this->covariance = $this->_dataCovariance * pow($this->factor, 2);
        $this->invCov = $this->_dataInvCov / pow($this->factor, 2);
        $this->normFactor = sqrt(2 * M_PI * $this->covariance) * count($this->dataset);
    }
214975
public void computeErrorVector( DMatrixRMaj A , DMatrixRMaj errors ) {
		errors.reshape(A.numRows,1);
		CommonOps_DDRM.mult(A,vectorT,errors);
	}
230977
def build_regression_matrix(H, model, build=None):
    
    ListOfTokens = model.split(' ')
    if H.shape[1]==1:
        size_index = len(str(H.shape[0]))
    else:
        size_index = len(str(H.shape[1]))
    
    if build is None:
        build = [True]*len(ListOfTokens)
    
    # Test if the vector has the wrong direction (lines instead of columns)
    if H.shape[0]==1:
        H = H.T
    
    # Collect the list of monomials
    Monom_Index = []
    for i in range(len(ListOfTokens)):
        if build[i]:
            Monom_Index += [grep(ListOfTokens, 'x' + str(0)*(size_index - \
            len(str(i))) + str(i))]
    
    Monom_Index = -np.sort(-Monom_Index)
    Monom_Index = np.unique(Monom_Index)
    
    if H.shape[1]==1:
        nb_var = H.shape[0]  # vector "mode": the number of vars is equal to the number of lines of H
        VectorMode = True
        
        for i in range(nb_var):
            for j in range(ListOfTokens.shape[0]):
                ListOfTokens[j] = ListOfTokens[j].replace(
                    'x' + str(0)*(size_index - len(str(i))) + str(i),
                    'H(' + str(i) + ')')
    else:
        nb_var = H.shape[0]  # matrix "mode": the number of vars is equal to the number of columns of H
        VectorMode = False
        
        for i in range(nb_var):
            for j in range(ListOfTokens.shape[0]):
                ListOfTokens[j] = ListOfTokens[j].replace(
                    'x' + str(0)*(size_index - len(str(i))) + str(i),
                    'H[i,' + str(i) + ')')
    
    # Now build the regression matrix
    if VectorMode:
        R = np.zeros((len(ListOfTokens), 1))
        for j in range(len(ListOfTokens)):
            R[j, 0] = eval(ListOfTokens[j])
    else:
        R = np.zeros((H.shape[0], len(ListOfTokens)))
        for i in range(H.shape[0]):
            for j in range(len(ListOfTokens)):
                R[i, j] = eval(ListOfTokens[j])
    
    return R
578443
def error_perturbation_var(C, S):
    

    """
    K = S.shape[0]
    cov = tmatrix_cov(C)
    for i in range(K):
        R = S[i, :, :]
        X[i] = (R[:, :, np.newaxis] * cov * R[:, np.newaxis, :]).sum()
225334
function() {
        var tmp = [], //12
			src = [], //16
			dst = [], //16
			mat = this.toArray(),
			i = 0;

        for (i = 0; i < 4; i++) {
            var i4 = i << 2;
            src[i] = mat[i4];
            src[i + 4] = mat[i4 + 1];
            src[i + 8] = mat[i4 + 2];
            src[i + 12] = mat[i4 + 3];
        }

        // calculate pairs for first 8 elements (cofactors)
        tmp[0] = src[10] * src[15];
        tmp[1] = src[11] * src[14];
        tmp[2] = src[9] * src[15];
        tmp[3] = src[11] * src[13];
        tmp[4] = src[9] * src[14];
        tmp[5] = src[10] * src[13];
        tmp[6] = src[8] * src[15];
        tmp[7] = src[11] * src[12];
        tmp[8] = src[8] * src[14];
        tmp[9] = src[10] * src[12];
        tmp[10] = src[8] * src[13];
        tmp[11] = src[9] * src[12];

        // calculate first 8 elements (cofactors)
        var src0 = src[0],
			src1 = src[1],
			src2 = src[2],
			src3 = src[3],
			src4 = src[4],
			src5 = src[5],
			src6 = src[6],
			src7 = src[7];
		dst[0] = tmp[0] * src5 + tmp[3] * src6 + tmp[4] * src7;
		dst[0] -= tmp[1] * src5 + tmp[2] * src6 + tmp[5] * src7;
		dst[1] = tmp[1] * src4 + tmp[6] * src6 + tmp[9] * src7;
		dst[1] -= tmp[0] * src4 + tmp[7] * src6 + tmp[8] * src7;
		dst[2] = tmp[2] * src4 + tmp[7] * src5 + tmp[10] * src7;
		dst[2] -= tmp[3] * src4 + tmp[6] * src5 + tmp[11] * src7;
		dst[3] = tmp[5] * src4 + tmp[8] * src5 + tmp[11] * src6;
		dst[3] -= tmp[4] * src4 + tmp[9] * src5 + tmp[10] * src6;
		dst[4] = tmp[1] * src1 + tmp[2] * src2 + tmp[5] * src3;
		dst[4] -= tmp[0] * src1 + tmp[3] * src2 + tmp[4] * src3;
		dst[5] = tmp[0] * src0 + tmp[7] * src2 + tmp[8] * src3;
		dst[5] -= tmp[1] * src0 + tmp[6] * src2 + tmp[9] * src3;
		dst[6] = tmp[3] * src0 + tmp[6] * src1 + tmp[11] * src3;
		dst[6] -= tmp[2] * src0 + tmp[7] * src1 + tmp[10] * src3;
		dst[7] = tmp[4] * src0 + tmp[9] * src1 + tmp[10] * src2;
		dst[7] -= tmp[5] * src0 + tmp[8] * src1 + tmp[11] * src2;

        // calculate pairs for second 8 elements (cofactors)
		tmp[0] = src2 * src7;
		tmp[1] = src3 * src6;
		tmp[2] = src1 * src7;
		tmp[3] = src3 * src5;
		tmp[4] = src1 * src6;
		tmp[5] = src2 * src5;
		tmp[6] = src0 * src7;
		tmp[7] = src3 * src4;
		tmp[8] = src0 * src6;
		tmp[9] = src2 * src4;
		tmp[10] = src0 * src5;
		tmp[11] = src1 * src4;

        // calculate second 8 elements (cofactors)
		src0 = src[8];
		src1 = src[9];
		src2 = src[10];
		src3 = src[11];
		src4 = src[12];
		src5 = src[13];
		src6 = src[14];
		src7 = src[15];
		dst[8] = tmp[0] * src5 + tmp[3] * src6 + tmp[4] * src7;
		dst[8] -= tmp[1] * src5 + tmp[2] * src6 + tmp[5] * src7;
		dst[9] = tmp[1] * src4 + tmp[6] * src6 + tmp[9] * src7;
		dst[9] -= tmp[0] * src4 + tmp[7] * src6 + tmp[8] * src7;
		dst[10] = tmp[2] * src4 + tmp[7] * src5 + tmp[10] * src7;
		dst[10] -= tmp[3] * src4 + tmp[6] * src5 + tmp[11] * src7;
		dst[11] = tmp[5] * src4 + tmp[8] * src5 + tmp[11] * src6;
		dst[11] -= tmp[4] * src4 + tmp[9] * src5 + tmp[10] * src6;
		dst[12] = tmp[2] * src2 + tmp[5] * src3 + tmp[1] * src1;
		dst[12] -= tmp[4] * src3 + tmp[0] * src1 + tmp[3] * src2;
		dst[13] = tmp[8] * src3 + tmp[0] * src0 + tmp[7] * src2;
		dst[13] -= tmp[6] * src2 + tmp[9] * src3 + tmp[1] * src0;
		dst[14] = tmp[6] * src1 + tmp[11] * src3 + tmp[3] * src0;
		dst[14] -= tmp[10] * src3 + tmp[2] * src0 + tmp[7] * src1;
		dst[15] = tmp[10] * src2 + tmp[4] * src0 + tmp[9] * src1;
		dst[15] -= tmp[8] * src1 + tmp[11] * src2 + tmp[5] * src0;

		var det = 1.0 / (src[0] * dst[0] + src[1] * dst[1] + src[2] * dst[2] + src[3] * dst[3]);
		for (i = 0, k = 0; i < 4; i++) {
			var m = this.matrix[i];
			for (var j = 0; j < 4; j++) {
				m[j] = dst[k++] * det;
			}
		}
		return this;
    }
1316296
def _get_covariance(self, X):
        
        result = pd.DataFrame(index=range(len(X)))
        column_names = self.get_column_names(X)
        for column_name in column_names:
            column = self.get_column(X, column_name)
            distrib = self.distribs[column_name]

            # get original distrib's cdf of the column
            cdf = distrib.cumulative_distribution(column)

            if distrib.constant_value is not None:
                # This is to avoid np.inf in the case the column is constant.
                cdf = np.ones(column.shape) - EPSILON

            # get inverse cdf using standard normal
            result = self.set_column(result, column_name, stats.norm.ppf(cdf))

        # remove any rows that have infinite values
        result = result[(result != np.inf).all(axis=1)]
        return pd.DataFrame(data=result).cov().values