docs: def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
          """
          See :meth:`superclass method
          <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
          for spec of input and result values.
          """
          C = self.COEFFS[imt]
  
          M = rup.mag - 6
          R = np.sqrt(dists.rjb ** 2 + C['h'] ** 2)
  
          # In the original formulation of the GMPE, distinction is only made
          # between rock and soil sites, which I assumed separated by the Vs30
          # value of 910m/s (see equation 5 of the paper)
          gamma = np.array([0 if v > 910. else 1 for v in sites.vs30])
  
          mean = np.zeros_like(R)
  
          mean += C['b1'] + \
                  C['b2'] * M + \
                  C['b3'] * M ** 2 + \
                  C['b5'] * np.log10(R) + \
                  C['b6'] * gamma
  
          # Convert from base 10 to base e
          mean /= np.log10(np.e)
  
          # Converting PSV to PSA
          if imt != PGA() and imt != PGV():
              omega = 2.*np.pi/imt.period
              mean += np.log(omega/(gravity*100))
  
          # Computing standard deviation
          stddevs = self._get_stddevs(C, stddev_types,  dists.rjb.shape[0])
  
          # Convert from base 10 to base e
          stddevs = [sd/np.log10(np.e) for sd in stddevs]
  
          return mean, stddevs