Will extract the standard deviations (volatilities) from the input covariance matrix. If "cleaning" is
enabled small variances will be replaced with a new minimal value.
515230
public static PrimitiveMatrix toVolatilities(Access2D<?> covariances, boolean clean) {

        int size = Math.toIntExact(Math.min(covariances.countRows(), covariances.countColumns()));

        PrimitiveMatrix.DenseReceiver retVal = PrimitiveMatrix.FACTORY.makeDense(size);

        if (clean) {

            MatrixStore<Double> covarianceMtrx = MatrixStore.PRIMITIVE.makeWrapper(covariances).get();

            double largest = covarianceMtrx.aggregateDiagonal(Aggregator.LARGEST);
            double limit = largest * size * PrimitiveMath.RELATIVELY_SMALL;
            double smallest = PrimitiveMath.SQRT.invoke(limit);

            for (int ij = 0; ij < size; ij++) {
                double variance = covariances.doubleValue(ij, ij);

                if (variance < limit) {
                    retVal.set(ij, smallest);
                } else {
                    retVal.set(ij, PrimitiveMath.SQRT.invoke(variance));
                }
            }

        } else {

            for (int ij = 0; ij < size; ij++) {
                double variance = covariances.doubleValue(ij, ij);

                if (variance <= PrimitiveMath.ZERO) {
                    retVal.set(ij, PrimitiveMath.ZERO);
                } else {
                    retVal.set(ij, PrimitiveMath.SQRT.invoke(variance));
                }
            }
        }

        return retVal.get();
    }
1295057
def resample(self):
        

        if self.just_resampled:
            warnings.warn(
                "Resampling without additional data; this may not perform as "
                "desired.",
                ResamplerWarning
            )

        # Record that we have performed a resampling step.
        self._just_resampled = True
        self._resample_count += 1

        # If we're tracking divergences, make a copy of the weights and
        # locations.
        if self._resampling_divergences is not None:
            old_locs = self.particle_locations.copy()
            old_weights = self.particle_weights.copy()

        # Record the previous mean, cov if needed.
        if self._debug_resampling:
            old_mean = self.est_mean()
            old_cov = self.est_covariance_mtx()

        # Find the new particle locations according to the chosen resampling
        # algorithm.
        # We pass the model so that the resampler can check for validity of
        # newly placed particles.
        # FIXME This feels fishy. If we update particles elsewwhere
        new_distribution = self.resampler(self.model, self)
        self.particle_weights = new_distribution.particle_weights
        self.particle_locations = new_distribution.particle_locations

        # Possibly canonicalize, if we've been asked to do so.
        if self._canonicalize:
            self.particle_locations[:, :] = self.model.canonicalize(self.particle_locations)

        # Instruct the model to clear its cache, demoting any errors to
        # warnings.
        try:
            self.model.clear_cache()
        except Exception as e:
            warnings.warn("Exception raised when clearing model cache: {}. Ignoring.".format(e))

        # Possibly track the new divergence.
        if self._resampling_divergences is not None:
            self._resampling_divergences.append(
                self._kl_divergence(old_locs, old_weights)
            )

        # Report current and previous mean, cov.
        if self._debug_resampling:
            new_mean = self.est_mean()
            new_cov = self.est_covariance_mtx()
            logger.debug("Resampling changed mean by {}. Norm change in cov: {}.".format(
                old_mean - new_mean,
                np.linalg.norm(new_cov - old_cov)
            ))
999471
public static double[] multinomialGaussianSample(double[] mean, double[][] covariance) {
        MultivariateNormalDistribution gaussian = new 
            MultivariateNormalDistribution(mean, covariance);
        gaussian.reseedRandomGenerator(RandomGenerator.getThreadLocalRandom().nextLong());
        return gaussian.sample();
    }
469422
def _get_stddevs(self, coeffs, stddev_types, num_sites):
        
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            stddevs.append(coeffs['sigma'] + np.zeros(num_sites))
        return np.array(stddevs)
1242754
def savgol_filter(x, window_length, polyorder, deriv=0, delta=1.0, axis=-1, mode='interp', cval=0.0):
    '''
    
    '''
    # linearly interpolate missing values before filtering
    x = np.ma.masked_invalid(pd.Series(x).interpolate())

    try:
        # start filtering from the first non-zero value since these won't be addressed by
        # the interpolation above
        ind = np.isfinite(x).nonzero()[0][0]
        x[ind:] = signal.savgol_filter(x[ind:], window_length, polyorder, deriv,
                                       delta, axis, mode, cval)
    except IndexError:
        pass
    return np.ma.masked_invalid(x)
388174
def rvs(self, size=1):
        
        return np.random.multivariate_normal(self.mean, self.cov, size)
47765
def shrink(self):
        
        x = np.nan_to_num(self.X.values)

        # de-mean returns
        t, n = np.shape(x)
        meanx = x.mean(axis=0)
        x = x - np.tile(meanx, (t, 1))

        # compute sample covariance matrix
        sample = (1.0 / t) * np.dot(x.T, x)

        # compute prior
        var = np.diag(sample).reshape(-1, 1)
        sqrtvar = np.sqrt(var)
        _var = np.tile(var, (n,))
        _sqrtvar = np.tile(sqrtvar, (n,))
        r_bar = (np.sum(sample / (_sqrtvar * _sqrtvar.T)) - n) / (n * (n - 1))
        prior = r_bar * (_sqrtvar * _sqrtvar.T)
        prior[np.eye(n) == 1] = var.reshape(-1)

        # compute shrinkage parameters and constant
        if self.delta is None:

            # what we call pi-hat
            y = x ** 2.0
            phi_mat = np.dot(y.T, y) / t - 2 * np.dot(x.T, x) * sample / t + sample ** 2
            phi = np.sum(phi_mat)

            # what we call rho-hat
            term1 = np.dot((x ** 3).T, x) / t
            help_ = np.dot(x.T, x) / t
            help_diag = np.diag(help_)
            term2 = np.tile(help_diag, (n, 1)).T * sample
            term3 = help_ * _var
            term4 = _var * sample
            theta_mat = term1 - term2 - term3 + term4
            theta_mat[np.eye(n) == 1] = np.zeros(n)
            rho = sum(np.diag(phi_mat)) + r_bar * np.sum(
                np.dot((1.0 / sqrtvar), sqrtvar.T) * theta_mat
            )

            # what we call gamma-hat
            gamma = np.linalg.norm(sample - prior, "fro") ** 2

            # compute shrinkage constant
            kappa = (phi - rho) / gamma
            shrinkage = max(0.0, min(1.0, kappa / t))
            self.delta = shrinkage
        else:
            # use specified constant
            shrinkage = self.delta

        # compute the estimator
        sigma = shrinkage * prior + (1 - shrinkage) * sample
        return self.format_and_annualise(sigma)
1113723
def sigmas_samples(self):
         """
        res = np.empty((self.nsamples, self.nstates, self.dimension), dtype=config.dtype)
        for i in range(self.nsamples):
            for j in range(self.nstates):
                res[i, j, :] = self._sampled_hmms[i].sigmas[j]
        return res
716979
def _covar_mstep_full(gmm, X, responsibilities, weighted_X_sum, norm,
                      min_covar):
    """"""
    # Eq. 12 from K. Murphy, "Fitting a Conditional Linear Gaussian
    # Distribution"
    n_features = X.shape[1]
    cv = np.empty((gmm.n_components, n_features, n_features))
    for c in range(gmm.n_components):
        post = responsibilities[:, c]
        mu = gmm.means_[c]
        diff = X - mu
        with np.errstate(under='ignore'):
            # Underflow Errors in doing post * X.T are  not important
            avg_cv = np.dot(post * diff.T, diff) / (post.sum() + 10 * EPS)
        cv[c] = avg_cv + min_covar * np.eye(n_features)
    return cv
938476
def get_deviations(self):
        
        bt = False
        if not self.istransformed:
            bt = True
            self._transform()
        mean_vec = self.mean()

        df = self.loc[:,:].copy()
        for col in df.columns:
            df.loc[:,col] -= mean_vec[col]
        if bt:
            self._back_transform()
        return type(self).from_dataframe(pst=self.pst,df=df)
269460
def shrink(self):
        
        x = np.nan_to_num(self.X.values)

        # de-mean returns
        t, n = np.shape(x)
        meanx = x.mean(axis=0)
        x = x - np.tile(meanx, (t, 1))
        xmkt = x.mean(axis=1).reshape(t, 1)

        # compute sample covariance matrix
        sample = np.cov(np.append(x, xmkt, axis=1), rowvar=False) * (t - 1) / t
        covmkt = sample[0:n, n].reshape(n, 1)
        varmkt = sample[n, n]
        sample = sample[:n, :n]
        prior = np.dot(covmkt, covmkt.T) / varmkt
        prior[np.eye(n) == 1] = np.diag(sample)

        # compute shrinkage parameters
        if self.delta is None:
            c = np.linalg.norm(sample - prior, "fro") ** 2
            y = x ** 2
            p = 1 / t * np.sum(np.dot(y.T, y)) - np.sum(sample ** 2)
            # r is divided into diagonal
            # and off-diagonal terms, and the off-diagonal term
            # is itself divided into smaller terms
            rdiag = 1 / t * np.sum(y ** 2) - sum(np.diag(sample) ** 2)
            z = x * np.tile(xmkt, (n,))
            v1 = 1 / t * np.dot(y.T, z) - np.tile(covmkt, (n,)) * sample
            roff1 = (
                np.sum(v1 * np.tile(covmkt, (n,)).T) / varmkt
                - np.sum(np.diag(v1) * covmkt.T) / varmkt
            )
            v3 = 1 / t * np.dot(z.T, z) - varmkt * sample
            roff3 = (
                np.sum(v3 * np.dot(covmkt, covmkt.T)) / varmkt ** 2
                - np.sum(np.diag(v3).reshape(-1, 1) * covmkt ** 2) / varmkt ** 2
            )
            roff = 2 * roff1 - roff3
            r = rdiag + roff

            # compute shrinkage constant
            k = (p - r) / c
            shrinkage = max(0, min(1, k / t))
            self.delta = shrinkage
        else:
            # use specified constant
            shrinkage = self.delta

        # compute the estimator
        sigma = shrinkage * prior + (1 - shrinkage) * sample
        return self.format_and_annualise(sigma)
1240193
def _extract_local_mean_gauss(image, mask = slice(None), sigma = 1, voxelspacing = None):
    
    # set voxel spacing
    if voxelspacing is None:
        voxelspacing = [1.] * image.ndim
        
    # determine gaussian kernel size in voxel units
    sigma = _create_structure_array(sigma, voxelspacing)
        
    return _extract_intensities(gaussian_filter(image, sigma), mask)
1165652
def _log_multivariate_normal_density_full(X, means, covars, min_covar=1.e-7):
    """"""
    n_samples, n_dim = X.shape
    nmix = len(means)
    log_prob = np.empty((n_samples, nmix))
    for c, (mu, cv) in enumerate(zip(means, covars)):
        try:
            cv_chol = linalg.cholesky(cv, lower=True)
        except linalg.LinAlgError:
            # The model is most probably stuck in a component with too
            # few observations, we need to reinitialize this components
            try:
                cv_chol = linalg.cholesky(cv + min_covar * np.eye(n_dim),
                                          lower=True)
            except linalg.LinAlgError:
                raise ValueError("'covars' must be symmetric, "
                                 "positive-definite")

        cv_log_det = 2 * np.sum(np.log(np.diagonal(cv_chol)))
        cv_sol = linalg.solve_triangular(cv_chol, (X - mu).T, lower=True).T
        log_prob[:, c] = - .5 * (np.sum(cv_sol ** 2, axis=1) +
                                 n_dim * np.log(2 * np.pi) + cv_log_det)

    return log_prob
442685
def predict_covariance(self, X, with_noise=True):
        
        _, v = self._predict(X, True, with_noise)
        return v
1241589
def min_volatility(self):
        """"""
        if not self.w:
            self.solve()
        var = []
        for w in self.w:
            a = np.dot(np.dot(w.T, self.cov_matrix), w)
            var.append(a)
        # return min(var)**.5, self.w[var.index(min(var))]
        self.weights = self.w[var.index(min(var))].reshape((self.n_assets,))
        return dict(zip(self.tickers, self.weights))
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)]
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)