This method proposes values for stochastics based on the empirical
        covariance of the values sampled so far.

        The proposal jumps are drawn from a multivariate normal distribution.
501716
def propose(self):
        

        arrayjump = np.dot(
            self.proposal_sd,
            np.random.normal(
                size=self.proposal_sd.shape[
                    0]))
        if self.verbose > 2:
            print_('Jump :', arrayjump)

        # Update each stochastic individually.
        for stochastic in self.stochastics:
            jump = arrayjump[self._slices[stochastic]].squeeze()
            if np.iterable(stochastic.value):
                jump = np.reshape(
                    arrayjump[
                        self._slices[
                            stochastic]],
                    np.shape(
                        stochastic.value))
            if self.isdiscrete[stochastic]:
                jump = round_array(jump)
            stochastic.value = stochastic.value + jump
1313529
def _get_warped_variance(self, mean, std, pred_init=None, deg_gauss_hermite=20):
        
        gh_samples, gh_weights = np.polynomial.hermite.hermgauss(deg_gauss_hermite)
        gh_samples = gh_samples[:, None]
        gh_weights = gh_weights[None, :]
        arg1 = gh_weights.dot(self._get_warped_term(mean, std, gh_samples, 
                                                    pred_init=pred_init) ** 2) / np.sqrt(np.pi)
        arg2 = self._get_warped_mean(mean, std, pred_init=pred_init,
                                     deg_gauss_hermite=deg_gauss_hermite)
        return arg1 - (arg2 ** 2)
518263
def _compute_stddevs(self, C, mag, rjb, imt, stddev_types):
        
        # aleatory uncertainty
        sigma_ale_m = np.interp(mag, [5.0, 5.5, 8.0],
                                [C['m50'], C['m55'], C['m80']])
        sigma_ale_rjb = np.interp(rjb, [5.0, 20.0], [C['r5'], C['r20']])
        sigma_ale = np.sqrt(sigma_ale_m ** 2 + sigma_ale_rjb ** 2)

        # epistemic uncertainty
        if imt.period < 1:
            sigma_epi = 0.36 + 0.07 * (mag - 6)
        else:
            sigma_epi = 0.34 + 0.06 * (mag - 6)

        sigma_total = np.sqrt(sigma_ale ** 2 + sigma_epi ** 2)

        stddevs = []
        for _ in stddev_types:
            stddevs.append(sigma_total)

        return stddevs
282451
function() {
			// Get the x and y coordinates
			var xCoords = this.pluck('x');
			var yCoords = this.pluck('y');

			// Get the means for the x and y coordinates
			var meanX = xCoords.mean();
			var meanY = yCoords.mean();

			// Get the standard deviations for the x and y coordinates
			var stdDevX = xCoords.stdDev();
			var stdDevY = yCoords.stdDev();

			// Map each element to the difference of the element and the mean
			// divided by the standard deviation
			xCoords.map(function(num) {
				return (num - meanX) / stdDevX;
			});
			yCoords.map(function(num) {
				return (num - meanY) / stdDevY;
			});

			// Multiply each element in the x by the corresponding value in
			// the y
			var nums = this.clone().map(function(num, index) {
				return xCoords.get(index) * yCoords.get(index);
			});

			// r is the sum of xCoords over the number of points minus 1
			return nums.sum() / (nums.size() - 1);
		}
393534
def covariance2(xs, ys)
      raise "Length xs and ys must be equal" unless xs.length == ys.length

      ev_x, ev_y = mean(xs), mean(ys)
      xs.zip(ys)
        .map{|x,y| (x-ev_x) * (y-ev_y)}
        .inject(0) {|sum, x| sum += x} / xs.length
    end
1133423
private int sampleFullConditional(int m, int n) {

        // remove z_i from the count variables
        int topic = z[m][n];
        word_topic_matrix[documents[m][n]][topic]--;
        nd[m][topic]--;
        nwsum[topic]--;
        ndsum[m]--;

        // do multinomial sampling via cumulative method:
        float[] p = new float[K];
        for (int k = 0; k < K; k++) {
            p[k] = (word_topic_matrix[documents[m][n]][k] + beta) / (nwsum[k] + V * beta)
                * (nd[m][k] + alpha) / (ndsum[m] + K * alpha);
        }
        
        topic = drawFromProbability(p);

        // add newly estimated z_i to count variables
        word_topic_matrix[documents[m][n]][topic]++;
        nd[m][topic]++;
        nwsum[topic]++;
        ndsum[m]++;

        return topic;
    }
1358914
public void process() {

		computeCovarince();

		float eigenvalue = smallestEigenvalue();
		// eigenvalue is the variance, convert to standard deviation
		double stdev = Math.sqrt(eigenvalue);
//		System.out.println("stdev "+stdev+"  total "+norm.size()+"  mean "+meanX+"  "+meanY);

		// approximate the spread in by doing it along the x-axis.
		// Really should be along the smallest singular axis
		double angle0 = Math.atan2(1.0,sigmas*(meanX-stdev));
		double angle1 = Math.atan2(1.0,sigmas*(meanX+stdev));

		spread = Math.abs(angle1-angle0);
	}
1027933
def stdevp(data, data_mean=None, data_variance=None):
    
    data_variance = data_variance or variance(data, data_mean)
    return sqrt(data_variance)
95576
def sample_gtf(data, D, k, likelihood='gaussian', prior='laplace',
                           lambda_hyperparams=None, lam_walk_stdev=0.01, lam0=1.,
                           dp_hyperparameter=None, w_hyperparameters=None,
                           iterations=7000, burn=2000, thin=10,
                           robust=False, empirical=False,
                           verbose=False):
    ''''''
    Dk = get_delta(D, k)
    dk_rows, dk_rowbreaks, dk_cols, dk_vals = decompose_delta(Dk)

    if likelihood == 'gaussian':
        y, w = data
    elif likelihood == 'binomial':
        trials, successes = data
    elif likelihood == 'poisson':
        obs = data
    else:
        raise Exception('Unknown likelihood type: {0}'.format(likelihood))

    if prior == 'laplace':
        if lambda_hyperparams == None:
            lambda_hyperparams = (1., 1.)
    elif prior == 'laplacegamma':
        if lambda_hyperparams == None:
            lambda_hyperparams = (1., 1.)
        if dp_hyperparameter == None:
            dp_hyperparameter = 1.
    elif prior == 'doublepareto' or prior == 'doublepareto2':
        if lambda_hyperparams == None:
            lambda_hyperparams = (1.0, 1.0)
        if dp_hyperparameter == None:
            dp_hyperparameter = 0.1
    elif prior == 'cauchy':
        if lambda_hyperparams == None:
            lambda_hyperparams = (1.0, 1.0)
    else:
        raise Exception('Unknown prior type: {0}.'.format(prior))

    if robust and w_hyperparameters is None:
        w_hyperparameters = (1., 1.)

    # Run the Gibbs sampler
    sample_size = (iterations - burn) / thin
    beta_samples = np.zeros((sample_size, D.shape[1]), dtype='double')
    lam_samples = np.zeros(sample_size, dtype='double')

    if likelihood == 'gaussian':
        if prior == 'laplace':
            gflbayes_gaussian_laplace(len(y), y, w,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'laplacegamma':
            if robust:
                gflbayes_gaussian_laplace_gamma_robust(len(y), y, w,
                                          dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                          lambda_hyperparams[0], lambda_hyperparams[1],
                                          dp_hyperparameter,
                                          w_hyperparameters[0], w_hyperparameters[1],
                                          iterations, burn, thin,
                                          double_matrix_to_c_pointer(beta_samples), lam_samples)
            else:    
                gflbayes_gaussian_laplace_gamma(len(y), y, w,
                                          dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                          lambda_hyperparams[0], lambda_hyperparams[1],
                                          dp_hyperparameter,
                                          iterations, burn, thin,
                                          double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'doublepareto':
            gflbayes_gaussian_doublepareto(len(y), y, w,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      lam_walk_stdev, lam0, dp_hyperparameter,
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'doublepareto2':
            gflbayes_gaussian_doublepareto2(len(y), y, w,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      dp_hyperparameter,
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'cauchy':
            gflbayes_gaussian_cauchy(len(y), y, w,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      lam_walk_stdev, lam0,
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
    elif likelihood == 'binomial':
        if prior == 'laplace':
            gflbayes_binomial_laplace(len(trials), trials, successes,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'doublepareto':
            gflbayes_binomial_doublepareto(len(trials), trials, successes,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      lam_walk_stdev, lam0, dp_hyperparameter,
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'laplacegamma':
            if empirical:
                gflbayes_empirical_binomial_laplace_gamma(len(trials), trials, successes,
                                          dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                          lam0,
                                          iterations, burn, thin,
                                          double_matrix_to_c_pointer(beta_samples), lam_samples)
            else:
                gflbayes_binomial_laplace_gamma(len(trials), trials, successes,
                                          dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                          lambda_hyperparams[0], lambda_hyperparams[1],
                                          dp_hyperparameter,
                                          iterations, burn, thin,
                                          double_matrix_to_c_pointer(beta_samples), lam_samples)
    elif likelihood == 'poisson':
        if prior == 'laplace':
            gflbayes_poisson_laplace(len(obs), obs,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
        elif prior == 'doublepareto':
            gflbayes_poisson_doublepareto(len(obs), obs,
                                      dk_rows, dk_rowbreaks, dk_cols, dk_vals,
                                      lambda_hyperparams[0], lambda_hyperparams[1],
                                      lam_walk_stdev, lam0, dp_hyperparameter,
                                      iterations, burn, thin,
                                      double_matrix_to_c_pointer(beta_samples), lam_samples)
    else:
        raise Exception('Unknown likelihood type: {0}'.format(likelihood))

    return (beta_samples,lam_samples)
753850
def _random_stochastic_matrix(m, n, k=None, sparse=False, format='csr',
                              random_state=None):
    
    if k is None:
        k = n
    # m prob vectors of dimension k, shape (m, k)
    probvecs = probvec(m, k, random_state=random_state)

    if k == n:
        P = probvecs
        if sparse:
            return scipy.sparse.coo_matrix(P).asformat(format)
        else:
            return P

    # if k < n:
    rows = np.repeat(np.arange(m), k)
    cols = \
        sample_without_replacement(
            n, k, num_trials=m, random_state=random_state
        ).ravel()
    data = probvecs.ravel()

    if sparse:
        P = scipy.sparse.coo_matrix((data, (rows, cols)), shape=(m, n))
        return P.asformat(format)
    else:
        P = np.zeros((m, n))
        P[rows, cols] = data
        return P
711538
def covariance(xs, ys)
      raise "Length xs and ys must be equal" unless xs.length == ys.length

      xys = xs.zip(ys).map{|x,y| x * y }
      mean(xys) - mean(xs) * mean(ys)
    end
1081197
def compute_mu_star_confidence(ee, num_trajectories, num_resamples,
                               conf_level):
    '''
    
    '''
    ee_resampled = np.zeros([num_trajectories])
    mu_star_resampled = np.zeros([num_resamples])

    if not 0 < conf_level < 1:
        raise ValueError("Confidence level must be between 0-1.")

    resample_index = np.random.randint(
        len(ee), size=(num_resamples, num_trajectories))
    ee_resampled = ee[resample_index]
    # Compute average of the absolute values over each of the resamples
    mu_star_resampled = np.average(np.abs(ee_resampled), axis=1)

    return norm.ppf(0.5 + conf_level / 2) * mu_star_resampled.std(ddof=1)
42562
def multi_reciprocal_extra(xs, ys, noise=False):
    
    ns = np.linspace(0.5, 6.0, num=56)
    best = ['', np.inf]
    fit_results = {}
    weights = get_weights(xs, ys)
    for n in ns:
        popt = extrapolate_reciprocal(xs, ys, n, noise)
        m = measure(reciprocal, xs, ys, popt, weights)
        pcov = []
        fit_results.update({n: {'measure': m, 'popt': popt, 'pcov': pcov}})
    for n in fit_results:
        if fit_results[n]['measure'] <= best[1]:
            best = reciprocal, fit_results[n]['measure'], n
    return fit_results[best[2]]['popt'], fit_results[best[2]]['pcov'], best
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
355928
def _simulate_coef_from_bootstraps(
            self, n_draws, coef_bootstraps, cov_bootstraps):
        """"""
        # Sample indices uniformly from {0, ..., n_bootstraps - 1}
        # (Wood pg. 199 step 6)
        random_bootstrap_indices = np.random.choice(
            np.arange(len(coef_bootstraps)), size=n_draws, replace=True)

        # Simulate `n_draws` many random coefficient vectors from a
        # multivariate normal distribution with mean and covariance given by
        # the bootstrap samples (indexed by `random_bootstrap_indices`) of
        # `coef_bootstraps` and `cov_bootstraps`. Because it's faster to draw
        # many samples from a certain distribution all at once, we make a dict
        # mapping bootstrap indices to draw indices and use the `size`
        # parameter of `np.random.multivariate_normal` to sample the draws
        # needed from that bootstrap sample all at once.
        bootstrap_index_to_draw_indices = defaultdict(list)
        for draw_index, bootstrap_index in enumerate(random_bootstrap_indices):
            bootstrap_index_to_draw_indices[bootstrap_index].append(draw_index)

        coef_draws = np.empty((n_draws, len(self.coef_)))

        for bootstrap, draw_indices in bootstrap_index_to_draw_indices.items():
            coef_draws[draw_indices] = np.random.multivariate_normal(
                coef_bootstraps[bootstrap], cov_bootstraps[bootstrap],
                size=len(draw_indices))

        return coef_draws
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
1136996
function sample_correlation(x, y) {
        var cov = sample_covariance(x, y),
            xstd = sample_standard_deviation(x),
            ystd = sample_standard_deviation(y);

        if (cov === null || xstd === null || ystd === null) {
            return null;
        }

        return cov / xstd / ystd;
    }