Returns the sigma mu adjustment factor
354821
def get_sigma_mu_adjustment(self, C, imt, rup, dists):
        
        if imt.name in "PGA PGV":
            # PGA and PGV are 2D arrays of dimension [nmags, ndists]
            sigma_mu = getattr(self, imt.name.lower())
            if rup.mag <= self.mags[0]:
                sigma_mu_m = sigma_mu[0, :]
            elif rup.mag >= self.mags[-1]:
                sigma_mu_m = sigma_mu[-1, :]
            else:
                intpl1 = interp1d(self.mags, sigma_mu, axis=0)
                sigma_mu_m = intpl1(rup.mag)
            # Linear interpolation with distance
            intpl2 = interp1d(self.dists, sigma_mu_m, bounds_error=False,
                              fill_value=(sigma_mu_m[0], sigma_mu_m[-1]))
            return intpl2(dists.rjb)
        # In the case of SA the array is of dimension [nmags, ndists, nperiods]
        # Get values for given magnitude
        if rup.mag <= self.mags[0]:
            sigma_mu_m = self.s_a[0, :, :]
        elif rup.mag >= self.mags[-1]:
            sigma_mu_m = self.s_a[-1, :, :]
        else:
            intpl1 = interp1d(self.mags, self.s_a, axis=0)
            sigma_mu_m = intpl1(rup.mag)
        # Get values for period - N.B. ln T, linear sigma mu interpolation
        if imt.period <= self.periods[0]:
            sigma_mu_t = sigma_mu_m[:, 0]
        elif imt.period >= self.periods[-1]:
            sigma_mu_t = sigma_mu_m[:, -1]
        else:
            intpl2 = interp1d(np.log(self.periods), sigma_mu_m, axis=1)
            sigma_mu_t = intpl2(np.log(imt.period))
        intpl3 = interp1d(self.dists, sigma_mu_t, bounds_error=False,
                          fill_value=(sigma_mu_t[0], sigma_mu_t[-1]))
        return intpl3(dists.rjb)
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
1205698
private double[] findAndEvaluateThreshold(DoubleDynamicHistogram histogram) {
    int n = histogram.getNumBins();
    double[] p1 = new double[n];
    double[] p2 = new double[n];
    double[] mu1 = new double[n];
    double[] mu2 = new double[n];
    double[] sigma1 = new double[n];
    double[] sigma2 = new double[n];
    double[] jt = new double[n];
    // Forward pass
    {
      MeanVariance mv = new MeanVariance();
      DoubleHistogram.Iter forward = histogram.iter();
      for(int i = 0; forward.valid(); i++, forward.advance()) {
        p1[i] = forward.getValue() + ((i > 0) ? p1[i - 1] : 0);
        mv.put(i, forward.getValue());
        mu1[i] = mv.getMean();
        sigma1[i] = mv.getNaiveStddev();
      }
    }
    // Backwards pass
    {
      MeanVariance mv = new MeanVariance();
      DoubleHistogram.Iter backwards = histogram.iter();
      backwards.seek(histogram.getNumBins() - 1); // Seek to last

      for(int j = n - 1; backwards.valid(); j--, backwards.retract()) {
        p2[j] = backwards.getValue() + ((j + 1 < n) ? p2[j + 1] : 0);
        mv.put(j, backwards.getValue());
        mu2[j] = mv.getMean();
        sigma2[j] = mv.getNaiveStddev();
      }
    }

    for(int i = 0; i < n; i++) {
      jt[i] = 1.0 + 2 * (p1[i] * (FastMath.log(sigma1[i]) - FastMath.log(p1[i])) + p2[i] * (FastMath.log(sigma2[i]) - FastMath.log(p2[i])));
    }

    int bestpos = -1;
    double bestgoodness = Double.NEGATIVE_INFINITY;

    double devPrev = jt[1] - jt[0];
    for(int i = 1; i < jt.length - 1; i++) {
      double devCur = jt[i + 1] - jt[i];
      // Local minimum found - calculate depth
      if(devCur >= 0 && devPrev <= 0) {
        double lowestMaxima = Double.POSITIVE_INFINITY;
        for(int j = i - 1; j > 0; j--) {
          if(jt[j - 1] < jt[j]) {
            lowestMaxima = Math.min(lowestMaxima, jt[j]);
            break;
          }
        }
        for(int j = i + 1; j < n - 2; j++) {
          if(jt[j + 1] < jt[j]) {
            lowestMaxima = Math.min(lowestMaxima, jt[j]);
            break;
          }
        }
        double localDepth = lowestMaxima - jt[i];

        final double mud = mu1[i] - mu2[i];
        double discriminability = mud * mud / (sigma1[i] * sigma1[i] + sigma2[i] * sigma2[i]);
        if(Double.isNaN(discriminability)) {
          discriminability = -1;
        }
        double goodness = localDepth * discriminability;
        if(goodness > bestgoodness) {
          bestgoodness = goodness;
          bestpos = i;
        }
      }
      devPrev = devCur;
    }
    Iter iter = histogram.iter();
    iter.seek(bestpos);
    return new double[] { iter.getRight(), bestgoodness };
  }
864005
def _sigmainf(N, h, m, dW, Km0, Pm0):
    """"""
    M = m*(m-1)//2
    Im = broadcast_to(np.eye(m), (N, m, m))
    IM = broadcast_to(np.eye(M), (N, M, M))
    Ims0 = np.eye(m**2)
    factor1 = broadcast_to((2.0/h)*np.dot(Km0, Ims0 - Pm0), (N, M, m**2))
    factor2 = _kp2(Im, _dot(dW, _t(dW)))
    factor3 = broadcast_to(np.dot(Ims0 - Pm0, Km0.T), (N, m**2, M))
    return 2*IM + _dot(_dot(factor1, factor2), factor3)
989920
def Normal(mu, sigma, tag=None):
    
    assert sigma > 0, 'Normal "sigma" must be greater than zero'
    return uv(ss.norm(loc=mu, scale=sigma), tag=tag)
446669
def _get_stddevs(self, C, stddev_types, rup, imt, num_sites):
        
        stddevs = []
        for stddev_type in stddev_types:
            sigma_mean = self._compute_standard_dev(rup, imt, C)
            sigma_tot = np.sqrt((sigma_mean ** 2) + (C['SigmaReg'] ** 2))
            sigma_tot = np.log10(np.exp(sigma_tot))
            stddevs.append(sigma_tot + np.zeros(num_sites))
        return stddevs
252498
public static double pdf(double x, double mu, double sigma) {
    if(x <= 0.) {
      return 0.;
    }
    final double xrel = (FastMath.log(x) - mu) / sigma;
    return 1 / (MathUtil.SQRTTWOPI * sigma * x) * FastMath.exp(-.5 * xrel * xrel);
  }
1183718
def get_stddevs(self, mean, stddev_types):
        
        mean = np.exp(mean)
        sigma = 0.39 + np.zeros(mean.shape)
        sigma[mean < 0.068] = 0.55
        idx = np.logical_and(mean >= 0.068, mean <= 0.21)
        sigma[idx] = 0.173- 0.140 * np.log(mean[idx])
        stddevs = []
        for stddev in stddev_types:
            if stddev == const.StdDev.TOTAL:
                stddevs.append(sigma)
        return stddevs
144133
def sigma_s2(self, r, R, r_ani, a, gamma, rho0_r0_gamma):
        
        beta = self._beta_ani(r, r_ani)
        return (1 - beta * R**2/r**2) * self.sigma_r2(r, a, gamma, rho0_r0_gamma, r_ani)
1363963
def norm(x, mu, sigma=1.0):
    
    return stats.norm(loc=mu, scale=sigma).pdf(x)
1221958
def dim_dtau(self, pars):
        
        """
        self._set_parameters(pars)
        # term 1
        num1 = -self.m * (self.w ** self.c) * self.c\
            * (self.tau ** (self.c - 1)) * np.sin(self.ang)
        term1 = self.sigmai * num1 / self.denom

        # term 2
        num2a = -self.m * self.otc * np.sin(self.ang)
        num2b = 2 * (self.w ** 2.0) * self.c * (self.tau ** (self.c - 1)) *\
            np.cos(self.ang)
        num2c = 2 * self.c * (self.w ** (self.c * 2)) *\
            (self.tau ** (2 * self.c - 1))
        term2 = self.sigma0 * num2a * (num2b + num2c) / (self.denom ** 2)

        result = term1 + term2

        return result
1126258
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]

        mean = (self.get_magnitude_scaling(C, rup.mag) +
                self.get_distance_term(C, rup, dists.rjb, imt) +
                self.get_site_amplification(C, sites))
        # GMPE originally in cm/s/s - convert to g
        if imt.name in "PGA SA":
            mean -= np.log(100.0 * g)
        stddevs = self.get_stddevs(C, dists.rjb.shape, stddev_types, sites)
        if self.sigma_mu_epsilon:
            # Apply the epistemic uncertainty factor (sigma_mu) multiplied by
            # the number of standard deviations
            sigma_mu = self.get_sigma_mu_adjustment(C, imt, rup, dists)
            # Cap sigma_mu at 0.5 ln units
            sigma_mu[sigma_mu > 0.5] = 0.5
            # Sigma mu should not be less than the standard deviation of the
            # fault-to-fault variability
            sigma_mu[sigma_mu < C["tau_fault"]] = C["tau_fault"]
            mean += (self.sigma_mu_epsilon * sigma_mu)
        return mean, stddevs
131766
def _get_stddevs(self, C, stddev_types, num_sites, mag_conversion_sigma):
        
        assert all(stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
                   for stddev_type in stddev_types)

        sigma = np.zeros(num_sites) + C['sigma'] * np.log(10)
        sigma = np.sqrt(sigma ** 2 + (C['a'] ** 2) * (mag_conversion_sigma ** 2))
        stddevs = [sigma for _ in stddev_types]

        return stddevs
708487
def dim_dm(self, pars):
        
        """
        self._set_parameters(pars)
        num1 = self.otc * np.sin(self.ang)
        result = -self.sigmai * num1 / self.denom

        return result
293024
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]
        mean = (self._get_magnitude_scaling(C, rup.mag) +
                self._get_distance_scaling(C, dists, rup.mag) +
                self._get_site_term(C, sites.vs30))

        # Mean is returned in terms of m/s^2. Need to convert to g
        mean -= np.log(g)
        stddevs = self.get_stddevs(C, sites.vs30.shape, stddev_types)
        return mean + self.adjustment_factor, stddevs
155400
def getModelPosterior(self,min):
        
        Sigma = self.getLaplaceCovar(min)
        n_params = self.vd.getNumberScales()
        ModCompl = 0.5*n_params*SP.log(2*SP.pi)+0.5*SP.log(SP.linalg.det(Sigma))
        RV = min['LML']+ModCompl
        return RV
892031
protected static double calcP_i(double f, double mu, double sigma) {
    final double fmu = f - mu;
    return ONEBYSQRT2PI / sigma * FastMath.exp(fmu * fmu / (-2 * sigma * sigma));
  }