INCREMENTAL MOVING SAMPLE VARIANCE // 
FUNCTION: incrmvariance( W )
Returns a method to compute a moving sample variance incrementally.

@param {Number} W - window size
@returns {Function} method to compute a moving sample variance incrementally
198858
function incrmvariance( W ) {
	if ( !isPositiveInteger( W ) ) {
		throw new TypeError( 'incrmvariance()::invalid input argument. Window size must be a positive integer. Value: `' + W + '`.' );
	}
	var arr = new Array( W ),
		n = W - 1,
		M2 = 0,
		mu = 0,
		N = 0,
		i = -1,
		delta,
		tmp,
		d1,
		d2;
	/**
	* FUNCTION: incrmvariance( [value] )
	*	If a `value` is provided, updates and returns the updated sample variance. If no `value` is provided, returns the current sample variance.
	*
	* @param {Number} [value] - value used to update the moving sample variance
	* @returns {Number|Null} sample variance or null
	*/
	return function incrmvariance( x ) {
		if ( !arguments.length ) {
			if ( N === 0 ) {
				return null;
			}
			if ( N === 1 ) {
				return 0;
			}
			if ( N < W ) {
				return M2 / (N-1);
			}
			return M2 / n;
		}
		// Update the index for managing the circular buffer...
		i = (i+1) % W;

		// Fill up the initial window; else, update the existing window...
		if ( N < W ) {
			arr[ i ] = x;
			N += 1;
			delta = x - mu;
			mu += delta / N;
			M2 += delta * (x - mu);
			if ( N === 1 ) {
				return 0;
			}
			return M2 / (N-1);
		}
		if ( N === 1 ) {
			return 0;
		}
		tmp = arr[ i ];
		arr[ i ] = x;
		delta = x - tmp;
		d1 = tmp - mu;
		mu += delta / W;
		d2 = x - mu;
		M2 += delta * (d1 + d2);

		return M2 / n;
	};
}
138549
public static double variance( GrayF64 img , double mean ) {

		if( BoofConcurrency.USE_CONCURRENT ) {
			return ImplImageStatistics_MT.variance(img,mean);
		} else {
			return ImplImageStatistics.variance(img,mean);
		}
	}
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
436528
def update(self, es, function_values, **kwargs):
        
        # TODO: on the linear function, the two mirrored samples lead
        # to a sharp increase of condition of the covariance matrix.
        # They should not be used to update the covariance matrix,
        # if the step-size inreases quickly. This should be fine with
        # negative updates though.
        if not self.initialized:
            self.initialize(es.N, es.opts)
        if 1 < 3:
            # use the ranking difference of the mirrors for adaptation
            # damp = 5 should be fine
            z = np.where(es.fit.idx == 1)[0][0] - np.where(es.fit.idx == 0)[0][0]
            z /= es.popsize - 1  # z in [-1, 1]
        self.s = (1 - self.sp.c) * self.s + self.sp.c * np.sign(z) * np.abs(z)**self.sp.z_exponent
        if self.s > 0:
            es.sigma *= exp(self.s / self.sp.dampup)
        else:
            es.sigma *= exp(self.s / self.sp.dampdown)
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]);
            }           
            
        }
    }
588549
public static double xbarVariance(TransposeDataCollection sampleDataCollection, int populationM, double Nbar) {
        double xbarVariance = 0.0;
        
        int sampleM = sampleDataCollection.size();
        
        double mean = mean(sampleDataCollection);
        
        for(Map.Entry<Object, FlatDataCollection> entry : sampleDataCollection.entrySet()) {
            double sum = 0.0;
            Iterator<Double> it = entry.getValue().iteratorDouble();
            while(it.hasNext()) {
                sum+=(it.next() - mean);
            }
            xbarVariance+= sum*sum/(sampleM-1);
        }
        
        xbarVariance *= (populationM-sampleM)/(populationM*sampleM*Nbar*Nbar);
        
        return xbarVariance;
    }
1244993
public static <T extends Number> WindowOver<T> variance(Expression<T> expr) {
        return new WindowOver<T>(expr.getType(), SQLOps.VARIANCE, expr);
    }
894390
function variance(values) {
	let avg = mean(values);
	return mean(values.map((e) => Math.pow(e - avg, 2)));
}
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
1347281
public static double variance( GrayS64 img , double mean ) {

		if( BoofConcurrency.USE_CONCURRENT ) {
			return ImplImageStatistics_MT.variance(img,mean);
		} else {
			return ImplImageStatistics.variance(img,mean);
		}
	}
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;		
	}
}
957565
def variance(array):
    
    if len(array) < 2:
        return 0
    u = mean(array)
    return sum([(x - u) ** 2 for x in array]) / (len(array) - 1)
664022
def moving_average(self, sampling_period,
                       window_size=None,
                       start=None, end=None,
                       placement='center',
                       pandas=False):
        
        start, end, mask = self._check_boundaries(start, end)

        # default to sampling_period if not given
        if window_size is None:
            window_size = sampling_period

        sampling_period = \
            self._check_regularization(start, end, sampling_period)

        # convert to datetime if the times are datetimes
        full_window = window_size * 1.  # convert to float if int or do nothing
        half_window = full_window / 2.  # divide by 2
        if (isinstance(start, datetime.datetime) and
                not isinstance(full_window, datetime.timedelta)):
            half_window = datetime.timedelta(seconds=half_window)
            full_window = datetime.timedelta(seconds=full_window)

        result = []
        current_time = start
        while current_time <= end:

            if placement == 'center':
                window_start = current_time - half_window
                window_end = current_time + half_window
            elif placement == 'left':
                window_start = current_time
                window_end = current_time + full_window
            elif placement == 'right':
                window_start = current_time - full_window
                window_end = current_time
            else:
                msg = 'unknown placement "{}"'.format(placement)
                raise ValueError(msg)

            # calculate mean over window and add (t, v) tuple to list
            try:
                mean = self.mean(window_start, window_end)
            except TypeError as e:
                if 'NoneType' in str(e):
                    mean = None
                else:
                    raise e
            result.append((current_time, mean))

            current_time += sampling_period

        # convert to pandas Series if pandas=True
        if pandas:

            try:
                import pandas as pd
            except ImportError:
                msg = "can't have pandas=True if pandas is not installed"
                raise ImportError(msg)

            result = pd.Series(
                [v for t, v in result],
                index=[t for t, v in result],
            )

        return result
521183
def variance
      μ = mean
      reduce(0) { |m, i| m + (i - μ).square } / (length - 1)
    end
1025776
def explained_variance(self):
        
        v = N.diagonal(self.covariance_matrix)
        return v[0:2].sum()/v.sum()
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
1013285
public void update(double sample) {
        // Adjust min&max
        min = (sample < min ? sample : min);
        max = (sample > max ? sample : max);

        // Update average
        average += (sample - average) / ++count;
        _pwrSumAverage += ( sample * sample - _pwrSumAverage) / count;

        // Update variance
        stdDev = Math.sqrt((_pwrSumAverage * count - count * average * average) / (count - 1));
        cv = 100 * (stdDev / average);
    }