Covariance matrix of the time shifted data
760391
def Ctt_(self):
        
        self._check_estimated()
        return self._rc.cov_YY(bessel=self.bessel)
54260
def coordinate_dimension_matrix(nc):
    '''
    
    '''
    retval = {}
    x = get_lon_variable(nc)
    if x:
        retval['x'] = nc.variables[x].dimensions
    y = get_lat_variable(nc)
    if y:
        retval['y'] = nc.variables[y].dimensions

    z = get_z_variable(nc)
    if z:
        retval['z'] = nc.variables[z].dimensions

    t = get_time_variable(nc)
    if t:
        retval['t'] = nc.variables[t].dimensions
    return retval
354354
private double[][] createMatrix() {
        double time = 0;
        double tmin = 0;
        tmin = rainData[0][0];
        time = tmin;
        double[][] matrix = new double[nTime][networkPipes.length + 1];

        for( int i = 0; i < nTime; ++i ) {
            matrix[i][0] = time;
            time += dt;
        }
        return matrix;
    }
6915
def _get_fit(self, time: int) -> typing.Tuple[np.ndarray, np.ndarray, float, float]:
        
        rawdata = self.averagedata[:, time]
        domain = np.arange(len(rawdata))
        datalength = len(domain)
        coefficients = np.zeros((datalength, self.function_number + 2))
        coefficients[:, 0] = 1
        coefficients[:, 1] = domain
        for i in range(self.function_number):
            coefficients[:, 2 + i] = self._gaussian_function(datalength, domain, 1, i)
        betas = linalg.inv(coefficients.transpose().dot(coefficients)).dot(coefficients.transpose().dot(rawdata))
        predicted_values = coefficients.dot(betas)
        residuals = rawdata - predicted_values
        error = np.sqrt(residuals.transpose().dot(residuals) / (datalength - (self.function_number + 2)))
        return predicted_values, residuals, residuals.mean(), error
209790
def summarize(self, n_timescales_to_report=5):
        """"""
        nonzeros = np.sum(np.abs(self.eigenvectors_) > 0, axis=0)
        active = '[%s]' % ', '.join(['%d/%d' % (n, self.n_features) for n in nonzeros[:n_timescales_to_report]])

        return """K-sparse time-structure Independent Components Analysis (tICA)
------------------------------------------------------------------
n_components        : {n_components}
shrinkage           : {shrinkage}
lag_time            : {lag_time}
kinetic_mapping     : {kinetic_mapping}
n_features          : {n_features}
Top {n_timescales_to_report} timescales :
{timescales}
Top {n_timescales_to_report} eigenvalues :
{eigenvalues}
Number of active degrees of freedom:
{active}
""".format(n_components=self.n_components, shrinkage=self.shrinkage_, lag_time=self.lag_time,
           kinetic_mapping=self.kinetic_mapping,
           timescales=self.timescales_[:n_timescales_to_report], eigenvalues=self.eigenvalues_[:n_timescales_to_report],
           n_features=self.n_features, active=active, n_timescales_to_report=n_timescales_to_report)
1025776
def explained_variance(self):
        
        v = N.diagonal(self.covariance_matrix)
        return v[0:2].sum()/v.sum()
861195
def calculate_covariance_matrix(X):
    

    n_features = X.shape[1]    
    S = np.zeros((n_features, n_features))
    m = np.mean(X, axis=0).reshape(n_features, 1)

    for x in X:
        v = x.reshape(n_features, 1) - m
        S += v @ v.T

    return 1/(X.shape[0]-1) * S
675716
def get_autoregressive_bias(max_length: int, dtype: str = C.DTYPE_FP32) -> mx.sym.Symbol:
    
    length_array = mx.sym.arange(max_length, dtype=dtype)
    # matrix with lower triangle and main diagonal set to 0, upper triangle set to 1
    bias = mx.sym.broadcast_greater(mx.sym.reshape(length_array, shape=(1, -1)),
                                    mx.sym.reshape(length_array, shape=(-1, 1)))
    bias = bias * -C.LARGE_VALUES[dtype]
    bias = mx.sym.reshape(bias, shape=(1, max_length, max_length))
    return mx.sym.BlockGrad(bias)
167945
def sample_from_likelihood(self, n_timesteps=10):
        

        self.latent_state_sequences = lmap(
            lambda A: ltake(
                n_timesteps,
                iterate(
                    lambda s: pd.Series(A @ s.values, index=s.index), self.s0
                ),
            ),
            self.transition_matrix_collection,
        )

        self.observed_state_sequences = [
            [self.sample_observed_state(s) for s in latent_state_sequence]
            for latent_state_sequence in self.latent_state_sequences
        ]
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()
1097014
def Cinv(self):
        """"""
        try:
            return np.linalg.inv(self.c)
        except np.linalg.linalg.LinAlgError:
            print('Warning: non-invertible noise covariance matrix c.')
            return np.eye(self.c.shape[0])
395013
def skew_matrix(w):
    ''''''
    return np.array([[0, -w[2], w[1]],
                     [w[2], 0, -w[0]],
                     [-w[1], w[0], 0]])
237590
def analyze(self, time: int=None) -> typing.Tuple[float, float, int, int]:
        
        if not time:
            time = int(len(self.times) / 2)
        if self.domains[time] == 0:
            yhat, residuals, mean_residual, error = self._get_fit(time)
            yhat_p = self.ddiff(yhat)
            yhat_pp = self.ddiff(yhat_p)
            noise = self._get_noise(residuals)
            curvature = (1 / self.ratio) * (1 / len(yhat_pp)) * np.sqrt(si.simps(yhat_pp ** 2))
            rng = (self.ratio * (np.max(self.averagedata[:, time]) -
                                 np.min(self.averagedata[:, time])))
            dmn = self.ratio * len(self.averagedata[:, time])

            self.noises[time] = np.log10(noise)
            self.curves[time] = np.log10(curvature)
            self.ranges[time] = np.log10(rng)
            self.domains[time] = np.log10(dmn)
        return self.noises[time], self.curves[time], self.ranges[time], self.domains[time]
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
910968
function autocorrelation(signal) {
        
        var n=signal.length;
        
        var i;
        
        var n2=nextpow2(n);

	if (n===0){
	    n2=0;
	}
        
        var paddedSignal=[];
        
        if (n==n2){
            paddedSignal=signal;
        }
        else{
            for (i=0;i<n2;i++){
                if (i<n){
                    paddedSignal[i]=signal[i];
                }
                else{
                    paddedSignal[i]=0;
                }
            }
        }

        var acv=autocovariance(paddedSignal);
        
        var variance=acv[0];
        
        var acf=[];
        
        var i;
        
        if (variance>0){
            for (i=0; i<n; i++){
                acf[i]=acv[i]/variance;    
            }
        }
        else{
            for (i=0; i<n; i++){
                acf[i]=0;    
            }            
        }
        return acf;
    
        function autocovariance(signal) {
            
            var n=signal.length;
            
            if (n===0)
                return [];
            
            var mean=0;
            var i;
            for (i=0; i<n; i++){
                mean+=signal[i];
            }
            mean/=n;
                      
            var padded_signal=[];
            
            for (i=0; i<n; i++){
                padded_signal[i]=signal[i]-mean;
                padded_signal[n+i]=0;
            }
                        
            var ft_signal = fft(padded_signal);
                       
            var pseudo_powerSpectralDensity=[];
            
            for (i=0; i<ft_signal.length; i++){
                pseudo_powerSpectralDensity[i]=complexMultiply(ft_signal[i], complexConjugate(ft_signal[i]));
            }
           
            var pseudo_autocovariance = ifft(pseudo_powerSpectralDensity);
    
            var mask=[];
            
            for (i=0; i<n; i++){
                mask[i]=1;
                mask[i+n]=0;
            }
            
            var ft_mask = fft(mask);
    
            var mask_powerSpectralDensity=[];
            
            for (i=0; i<ft_mask.length; i++){
                mask_powerSpectralDensity[i]=complexMultiply(ft_mask[i], complexConjugate(ft_mask[i]));
            }
            
            var mask_correction_factors = ifft(mask_powerSpectralDensity);
    
            var acv=[];
            
            for (i=0; i<n; i++){
                acv[i]=complexDivideRealPart(pseudo_autocovariance[i],mask_correction_factors[i]);
            }    
            
            return acv;
            /*
             * ComplexMultiply
             * 
             * From complex.js in fft-js in https://github.com/vail-systems/node-fft
             *
             * (c) Vail Systems. Joshua Jung and Ben Bryan. 2015
             */
            function complexMultiply(a, b){
                return [(a[0] * b[0] - a[1] * b[1]), 
                        (a[0] * b[1] + a[1] * b[0])];
            }
            
            function complexConjugate(a){
                return [a[0], -a[1]];
            }
            /*function complexDivide(a,b){
                return [ (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]),
                         (a[1]*b[0]-a[0]*b[1])/(b[0]*b[0]+b[1]*b[1])];
            }*/
            //We only need the real part here
            function complexDivideRealPart(a,b){
                return (a[0]*b[0]+a[1]*b[1])/(b[0]*b[0]+b[1]*b[1]);
            }           
            
        }
    }
1102124
def util_pattern_space(time_series, lag, dim):
    
    n = len(time_series)

    if lag * dim > n:
        raise Exception('Result matrix exceeded size limit, try to change lag or dim.')
    elif lag < 1:
        raise Exception('Lag should be greater or equal to 1.')

    pattern_space = np.empty((n - lag * (dim - 1), dim))
    for i in range(n - lag * (dim - 1)):
        for j in range(dim):
            pattern_space[i][j] = time_series[i + j * lag]

    return pattern_space
569526
def corrgroups60(display=False):
    

    # set a constant seed
    old_seed = np.random.seed()
    np.random.seed(0)

    # generate dataset with known correlation
    N = 1000
    M = 60

    # set one coefficent from each group of 3 to 1
    beta = np.zeros(M)
    beta[0:30:3] = 1

    # build a correlation matrix with groups of 3 tightly correlated features
    C = np.eye(M)
    for i in range(0,30,3):
        C[i,i+1] = C[i+1,i] = 0.99
        C[i,i+2] = C[i+2,i] = 0.99
        C[i+1,i+2] = C[i+2,i+1] = 0.99
    f = lambda X: np.matmul(X, beta)

    # Make sure the sample correlation is a perfect match
    X_start = np.random.randn(N, M)
    X_centered = X_start - X_start.mean(0)
    Sigma = np.matmul(X_centered.T, X_centered) / X_centered.shape[0]
    W = np.linalg.cholesky(np.linalg.inv(Sigma)).T
    X_white = np.matmul(X_centered, W.T)
    assert np.linalg.norm(np.corrcoef(np.matmul(X_centered, W.T).T) - np.eye(M)) < 1e-6 # ensure this decorrelates the data

    # create the final data
    X_final = np.matmul(X_white, np.linalg.cholesky(C).T)
    X = X_final
    y = f(X) + np.random.randn(N) * 1e-2

    # restore the previous numpy random seed
    np.random.seed(old_seed)

    return pd.DataFrame(X), y