Returns the intra-event random effects coefficient (phi)
        Equation 28.
16027
def _get_philny(self, C, mag):
        
        if mag <= 4.5:
            return C["phi1"]
        elif mag >= 5.5:
            return C["phi2"]
        else:
            return C["phi2"] + (C["phi1"] - C["phi2"]) * (5.5 - mag)
297626
def _get_phi(self, C, mag):
        
        if mag < 5.5:
            return C["phi1"]
        elif mag < 5.75:
            return C["phi1"] + (C["phi2"] - C["phi1"]) * ((mag - 5.5) / 0.25)
        else:
            return C["phi2"]
630918
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        

        COEFFS = self.COEFFS[imt]
        R = self._compute_term_r(COEFFS, rup.mag, dists.rrup)

        mean = 10 ** (self._compute_mean(COEFFS, rup.mag, R))

        # Convert units to g,
        # but only for PGA and SA (not PGV):
        if imt.name in "SA PGA":
            mean = np.log(mean / (g*100.))
        else:
            # PGV:
            mean = np.log(mean)

        c1_rrup = _compute_C1_term(COEFFS, dists.rrup)
        log_phi_ss = 1.00

        stddevs = self._get_stddevs(
            COEFFS, stddev_types, sites.vs30.shape[0], rup.mag, c1_rrup,
            log_phi_ss, COEFFS['mean_phi_ss']
        )

        return mean, stddevs
909278
def _get_stddevs(self, C, C_PGA, rup, sites, pga1100, stddev_types):
        
        # Get stddevs for PGA on basement rock
        tau_lnpga_b, phi_lnpga_b = self._get_stddevs_pga(C_PGA, rup)
        num_sites = len(sites.vs30)
        # Get tau_lny on the basement rock
        tau_lnyb = self._get_taulny(C, rup.mag)
        # Get phi_lny on the basement rock
        phi_lnyb = np.sqrt(self._get_philny(C, rup.mag) ** 2. -
                           self.CONSTS["philnAF"] ** 2.)
        # Get site scaling term
        alpha = self._get_alpha(C, sites.vs30, pga1100)
        # Evaluate tau according to equation 29
        tau = np.sqrt(
            (tau_lnyb ** 2.) +
            ((alpha ** 2.) * (tau_lnpga_b ** 2.)) +
            (2.0 * alpha * C["rholny"] * tau_lnyb * tau_lnpga_b))

        # Evaluate phi according to equation 30
        phi = np.sqrt(
            (phi_lnyb ** 2.) +
            (self.CONSTS["philnAF"] ** 2.) +
            ((alpha ** 2.) * (phi_lnpga_b ** 2.)) +
            (2.0 * alpha * C["rholny"] * phi_lnyb * phi_lnpga_b))
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(np.sqrt((tau ** 2.) + (phi ** 2.)) +
                               np.zeros(num_sites))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(phi + np.zeros(num_sites))
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(tau + np.zeros(num_sites))
        return stddevs
386414
def _get_stddevs(self, C, rup, dists, sites, stddev_types):
        
        stddevs = []
        num_sites = len(sites.vs30)
        tau = self._get_inter_event_tau(C, rup.mag, num_sites)
        phi = self._get_intra_event_phi(C,
                                        rup.mag,
                                        dists.rjb,
                                        sites.vs30,
                                        num_sites)
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(np.sqrt((tau ** 2.0) + (phi ** 2.0)))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(phi)
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(tau)
        return stddevs
576185
def coefficient_of_variation(data, data_mean=None, data_stdev=None):
    
    data_mean = data_mean or mean(data)
    data_stdev = data_stdev or stdevp(data, data_mean)
    if data_mean == 0:
        return float("inf") if data_stdev != 0 else 0
    else:
        return data_stdev / data_mean
8486
def phi( n ):
  """"""

  assert isinstance( n, integer_types )

  if n < 3: return 1

  result = 1
  ff = factorization( n )
  for f in ff:
    e = f[1]
    if e > 1:
      result = result * f[0] ** (e-1) * ( f[0] - 1 )
    else:
      result = result * ( f[0] - 1 )
  return result
248809
def _get_stddevs(self, C, C_PGA, pga1100, rup, sites, stddev_types):
        
        std_intra = self._compute_intra_event_std(C, C_PGA, pga1100, rup.mag,
                                                  sites.vs30,
                                                  sites.vs30measured)
        std_inter = self._compute_inter_event_std(C, C_PGA, pga1100, rup.mag,
                                                  sites.vs30)
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(np.sqrt(std_intra ** 2 + std_inter ** 2))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(std_intra)
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(std_inter)
        return stddevs
407111
def _compute_inter_event_std(self, C, C_PGA, pga1100, mag, vs30):
        
        tau_0 = self._compute_std_0(C['s3'], C['s4'], mag)
        tau_b_pga = self._compute_std_0(C_PGA['s3'], C_PGA['s4'], mag)
        delta_amp = self._compute_partial_derivative_site_amp(C, pga1100, vs30)

        std_inter = np.sqrt(tau_0 ** 2 + (delta_amp ** 2) * (tau_b_pga ** 2) +
                            2 * delta_amp * tau_0 * tau_b_pga * C['rho'])

        return std_inter
1316296
def _get_covariance(self, X):
        
        result = pd.DataFrame(index=range(len(X)))
        column_names = self.get_column_names(X)
        for column_name in column_names:
            column = self.get_column(X, column_name)
            distrib = self.distribs[column_name]

            # get original distrib's cdf of the column
            cdf = distrib.cumulative_distribution(column)

            if distrib.constant_value is not None:
                # This is to avoid np.inf in the case the column is constant.
                cdf = np.ones(column.shape) - EPSILON

            # get inverse cdf using standard normal
            result = self.set_column(result, column_name, stats.norm.ppf(cdf))

        # remove any rows that have infinite values
        result = result[(result != np.inf).all(axis=1)]
        return pd.DataFrame(data=result).cov().values
1199316
def _get_SRF_phi(self, imt_per):
        
        if imt_per < 0.6:
            srf = 0.8
        elif 0.6 <= imt_per < 1:
            srf = self._interp_function(0.7, 0.8, 1, 0.6, imt_per)
        elif 1 <= imt_per <= 10:
            srf = self._interp_function(0.6, 0.7, 10, 1, imt_per)
        else:
            srf = 1

        return srf
1292308
def get_phi_ss_at_quantile(phi_model, quantile):
    
    # Setup SA coeffs - the backward compatible Python 2.7 way
    coeffs = deepcopy(phi_model.sa_coeffs)
    coeffs.update(phi_model.non_sa_coeffs)
    for imt in coeffs:
        if quantile is None:
            coeffs[imt] = {"a": phi_model[imt]["mean_a"],
                           "b": phi_model[imt]["mean_b"]}
        else:
            coeffs[imt] = {
                "a": _at_percentile(phi_model[imt]["mean_a"],
                                    phi_model[imt]["var_a"],
                                    quantile),
                "b": _at_percentile(phi_model[imt]["mean_b"],
                                    phi_model[imt]["var_b"],
                                    quantile)
                }
    return CoeffsTable(sa_damping=5., table=coeffs)
236214
def PC_varExplained(Y,standardized=True):
    
    # figuring out the number of latent factors
    if standardized:
        Y-=Y.mean(0)
        Y/=Y.std(0)
    covY = sp.cov(Y)
    S,U = linalg.eigh(covY+1e-6*sp.eye(covY.shape[0]))
    S = S[::-1]
    rv = np.array([S[0:i].sum() for i in range(1,S.shape[0])])
    rv/= S.sum()
    return rv
1246578
def _get_stddevs(self, C, imt, rup, sites, stddev_types, sa1180, dists):
        
        std_intra = self._get_intra_event_std(C, rup.mag, sa1180, sites.vs30,
                                              sites.vs30measured, dists.rrup)
        std_inter = self._get_inter_event_std(C, rup.mag, sa1180, sites.vs30)
        stddevs = []
        for stddev_type in stddev_types:
            assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
            if stddev_type == const.StdDev.TOTAL:
                stddevs.append(np.sqrt(std_intra ** 2 +
                                       std_inter ** 2))
            elif stddev_type == const.StdDev.INTRA_EVENT:
                stddevs.append(std_intra)
            elif stddev_type == const.StdDev.INTER_EVENT:
                stddevs.append(std_inter)
        return stddevs
142733
def moments_match_ep(self, Y_i, tau_i, v_i, Y_metadata_i=None):
        
        if Y_i == 1:
            sign = 1.
        elif Y_i == 0 or Y_i == -1:
            sign = -1
        else:
            raise ValueError("bad value for Bernoulli observation (0, 1)")
        if isinstance(self.gp_link, link_functions.Probit):
            z = sign*v_i/np.sqrt(tau_i**2 + tau_i)
            phi_div_Phi = derivLogCdfNormal(z)
            log_Z_hat = logCdfNormal(z)

            mu_hat = v_i/tau_i + sign*phi_div_Phi/np.sqrt(tau_i**2 + tau_i)
            sigma2_hat = 1./tau_i - (phi_div_Phi/(tau_i**2+tau_i))*(z+phi_div_Phi)

        elif isinstance(self.gp_link, link_functions.Heaviside):
            z = sign*v_i/np.sqrt(tau_i)
            phi_div_Phi = derivLogCdfNormal(z)
            log_Z_hat = logCdfNormal(z)
            mu_hat = v_i/tau_i + sign*phi_div_Phi/np.sqrt(tau_i)
            sigma2_hat = (1. - a*phi_div_Phi - np.square(phi_div_Phi))/tau_i
        else:
            #TODO: do we want to revert to numerical quadrature here?
            raise ValueError("Exact moment matching not available for link {}".format(self.gp_link.__name__))

        # TODO: Output log_Z_hat instead of Z_hat (needs to be change in all others likelihoods)
        return np.exp(log_Z_hat), mu_hat, sigma2_hat
1188283
def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
        
        # extracting dictionary of coefficients specific to required
        # intensity measure type.

        C = self.COEFFS[imt]
        imean = self._get_mean(C, rup, dists, sites)
        if imt.name in "SA PGA":
            # Convert units to g,
            # but only for PGA and SA (not PGV):
            mean = np.log((10.0 ** (imean - 2.0)) / g)
        else:
            # PGV:
            mean = np.log(10.0 ** imean)

        istddevs = self._get_stddevs(C, stddev_types, len(sites.vs30))
        stddevs = np.log(10.0 ** np.array(istddevs))
        return mean + self.adjustment_factor, stddevs
987207
@Override
	public RandomVariable getCovariance(int timeIndex, int component1, int component2, RandomVariable[] realizationAtTimeIndex) {

		RandomVariable volatilityOfComponent1 = volatilityModel.getVolatility(timeIndex, component1);
		RandomVariable volatilityOfComponent2 = volatilityModel.getVolatility(timeIndex, component2);

		double					correlationOfComponent1And2 = correlationModel.getCorrelation(timeIndex, component1, component2);

		RandomVariable covariance = volatilityOfComponent1.mult(volatilityOfComponent2).mult(correlationOfComponent1And2);

		return covariance;
	}