Compute the diagonal of the covariance matrix associated to X.
502133
def Kdiag(self, X):
        """"""
        vyt = self.variance_Yt
        vyx = self.variance_Yx

        lyt = 1./(2*self.lengthscale_Yt)
        lyx = 1./(2*self.lengthscale_Yx)

        a = self.a
        b = self.b
        c = self.c

        ## dk^2/dtdt'
        k1 = (2*lyt )*vyt*vyx
        ## dk^2/dx^2
        k2 = ( - 2*lyx )*vyt*vyx
        ## dk^4/dx^2dx'^2
        k3 = ( 4*3*lyx**2 )*vyt*vyx


        Kdiag = np.zeros(X.shape[0])
        slices = index_to_slices(X[:,-1])

        for i, ss1 in enumerate(slices):
            for s1 in ss1:
                if i==0:
                    Kdiag[s1]+= vyt*vyx
                elif i==1:
                    #i=1
                    Kdiag[s1]+= b**2*k1 - 2*a*c*k2 + a**2*k3 + c**2*vyt*vyx
                    #Kdiag[s1]+= Vu*Vy*(k1+k2+k3)
                else:
                    raise ValueError("invalid input/output index")

        return Kdiag
1300178
def diagonalize(self):
        ''''''
        self.eigvals, self.eigvecs = np.linalg.eig(
            (self.tensor.transpose() + self.tensor) / 2.0)
        self.eigvals = np.diag(np.dot(
            np.dot(self.eigvecs.transpose(), self.tensor), self.eigvecs))
352964
def r(w, Xs):
    '''
    '''
    n = len(Xs)
    P = projection_matrix(w)
    c = C(w, Xs)

    return np.sqrt(sum(np.dot(c - X, np.dot(P, c - X)) for X in Xs) / n)
1272726
def _covar_mstep_diag(gmm, X, responsibilities, weighted_X_sum, norm,
                      min_covar):
    """"""
    avg_X2 = np.dot(responsibilities.T, X * X) * norm
    avg_means2 = gmm.means_ ** 2
    avg_X_means = gmm.means_ * weighted_X_sum * norm
    return avg_X2 - 2 * avg_X_means + avg_means2 + min_covar
650799
def _distance_covariance_sqr_fast_generic(
        x, y, unbiased=False):  # pylint:disable=too-many-locals
    # This function has many locals so it can be compared
    # with the original algorithm.
    """"""
    x = np.asarray(x)
    y = np.asarray(y)

    x = np.ravel(x)
    y = np.ravel(y)

    n = x.shape[0]
    assert n > 3
    assert n == y.shape[0]
    temp = range(n)

    # Step 1
    ix0 = np.argsort(x)
    vx = x[ix0]

    ix = np.zeros(n, dtype=int)
    ix[ix0] = temp

    iy0 = np.argsort(y)
    vy = y[iy0]

    iy = np.zeros(n, dtype=int)
    iy[iy0] = temp

    # Step 2
    sx = np.cumsum(vx)
    sy = np.cumsum(vy)

    # Step 3
    alpha_x = ix
    alpha_y = iy

    beta_x = sx[ix] - vx[ix]
    beta_y = sy[iy] - vy[iy]

    # Step 4
    x_dot = np.sum(x)
    y_dot = np.sum(y)

    # Step 5
    a_i_dot = x_dot + (2 * alpha_x - n) * x - 2 * beta_x
    b_i_dot = y_dot + (2 * alpha_y - n) * y - 2 * beta_y

    sum_ab = np.sum(a_i_dot * b_i_dot)

    # Step 6
    a_dot_dot = 2 * np.sum(alpha_x * x) - 2 * np.sum(beta_x)
    b_dot_dot = 2 * np.sum(alpha_y * y) - 2 * np.sum(beta_y)

    # Step 7
    gamma_1 = _partial_sum_2d(x, y, np.ones(n, dtype=x.dtype))
    gamma_x = _partial_sum_2d(x, y, x)
    gamma_y = _partial_sum_2d(x, y, y)
    gamma_xy = _partial_sum_2d(x, y, x * y)

    # Step 8
    aijbij = np.sum(x * y * gamma_1 + gamma_xy - x * gamma_y - y * gamma_x)

    if unbiased:
        d3 = (n - 3)
        d2 = (n - 2)
        d1 = (n - 1)
    else:
        d3 = d2 = d1 = n

    # Step 9
    d_cov = (aijbij / n / d3 - 2 * sum_ab / n / d2 / d3 +
             a_dot_dot / n * b_dot_dot / d1 / d2 / d3)

    return d_cov
93397
def _calc_q_statistic(x, h, nt):
    
    t, m, n = x.shape

    # covariance matrix of x
    c0 = acm(x, 0)

    # LU factorization of covariance matrix
    c0f = sp.linalg.lu_factor(c0, overwrite_a=False, check_finite=True)

    q = np.zeros((3, h + 1))
    for l in range(1, h + 1):
        cl = acm(x, l)

        # calculate tr(cl' * c0^-1 * cl * c0^-1)
        a = sp.linalg.lu_solve(c0f, cl)
        b = sp.linalg.lu_solve(c0f, cl.T)
        tmp = a.dot(b).trace()

        # Box-Pierce
        q[0, l] = tmp

        # Ljung-Box
        q[1, l] = tmp / (nt - l)

        # Li-McLeod
        q[2, l] = tmp

    q *= nt
    q[1, :] *= (nt + 2)

    q = np.cumsum(q, axis=1)

    for l in range(1, h+1):
        q[2, l] = q[0, l] + m * m * l * (l + 1) / (2 * nt)

    return q
178673
def compute(self, x_arr, y_arr):
        '''
        
        '''
        return np.linalg.norm(x_arr - y_arr, axis=-1)
1153679
def cnst_A(self, X):
        
        """

        return np.concatenate(
            [sl.Gax(X, ax)[..., np.newaxis] for ax in self.axes] +
            [X[..., np.newaxis],], axis=X.ndim)
1058492
def dirichlet_covariance(alpha):
    

    """
    alpha0 = alpha.sum()
    norm = alpha0 ** 2 * (alpha0 + 1.0)

    """Non normalized covariance"""
    Z = -alpha[:, np.newaxis] * alpha[np.newaxis, :]

    """Correct diagonal"""
    ind = np.diag_indices(Z.shape[0])
    Z[ind] += alpha0 * alpha

    """Covariance matrix"""
    cov = Z / norm

    return cov
1025776
def explained_variance(self):
        
        v = N.diagonal(self.covariance_matrix)
        return v[0:2].sum()/v.sum()
381182
def gradients_X(self, dL_dK, X, X2, target):
        """"""
        # First account for gradients arising from presence of X in exponent.
        self._K_computations(X, X2)
        if X2 is None:
            _K_dist = 2*(X[:, None, :] - X[None, :, :])
        else:
            _K_dist = X[:, None, :] - X2[None, :, :] # don't cache this in _K_co
        gradients_X = (-2.*self.variance)*np.transpose((self._K_dvar/self._w2)[:, :, None]*_K_dist, (1, 0, 2))
        target += np.sum(gradients_X*dL_dK.T[:, :, None], 0)
        # Now account for gradients arising from presence of X in lengthscale.
        self._dK_computations(dL_dK)
        if X2 is None:
            target += 2.*self.mapping.df_dX(self._dL_dl[:, None], X)
        else:
            target += self.mapping.df_dX(self._dL_dl[:, None], X)
589311
def cnst_AT(self, X):
        
        """

        Xf = sl.rfftn(X, axes=self.axes)
        return np.sum(sl.irfftn(np.conj(self.GAf)*Xf, self.axsz,
                                axes=self.axes), axis=self.Y.ndim-1)
287606
def _pca(x, k=2):
    ""
    x = x-torch.mean(x,0)
    U,S,V = torch.svd(x.t())
    return torch.mm(x,U[:,:k])
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))
934353
def fit(self, X):
        
        U, V = self.split_matrix(X)
        self.tau = stats.kendalltau(U, V)[0]
        self.theta = self.compute_theta()
        self.check_theta()
741590
def Kdiag(self, X):
        """"""
        K_ff = np.ones(X.shape[0])
        for i,x in enumerate(X):
            for il,l in enumerate(self.lengthscale):
                idx = il*2
                K_ff[i] *= self.k_ff(x[idx],x[idx],l)
        return K_ff * self.variances[0]
231217
private static void gaussian(float a[][],
	                             int index[]) {
		int n = index.length;
		float c[] = new float[n];
		// Initialize the index
		for (int i = 0; i < n; ++i) {
			index[i] = i;
		}
		// Find the rescaling factors, one from each row
		for (int i = 0; i < n; ++i){
			float c1 = 0;
			for (int j = 0; j < n; ++j){
				float c0 = Math.abs(a[i][j]);
				if (c0 > c1) {
					c1 = c0;
				}
			}
			c[i] = c1;
		}
		// Search the pivoting element from each column
		int k = 0;
		for (int j = 0; j < n - 1; ++j){
			float pi1 = 0;
			for (int i = j; i < n; ++i){
				float pi0 = Math.abs(a[index[i]][j]);
				pi0 /= c[index[i]];
				if (pi0 > pi1) {

					pi1 = pi0;
					k = i;
				}
			}
			// Interchange rows according to the pivoting order
			int itmp = index[j];
			index[j] = index[k];
			index[k] = itmp;
			for (int i = j + 1; i < n; ++i){
				float pj = a[index[i]][j] / a[index[j]][j];
				// Record pivoting ratios below the diagonal
				a[index[i]][j] = pj;

				// Modify other elements accordingly
				for (int l = j + 1; l < n; ++l){
					a[index[i]][l] -= pj * a[index[j]][l];
				}
			}
		}
	}