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.
          """
          # extracting dictionary of coefficients specific to required
          # intensity measure type.
          C = self.COEFFS[imt]
          C_SR = self.COEFFS_SOIL_RESPONSE[imt]
  
          # compute PGA on rock conditions - needed to compute non-linear
          # site amplification term
          pga4nl = self._get_pga_on_rock(rup, dists, C)
  
          # equation 1, pag 106, without sigma term, that is only the first 3
          # terms. The third term (site amplification) is computed as given in
          # equation (6), that is the sum of a linear term - equation (7) - and
          # a non-linear one - equations (8a) to (8c).
          # Mref, Rref values are given in the caption to table 6, pag 119.
          if imt == PGA():
              # avoid recomputing PGA on rock, just add site terms
              mean = np.log(pga4nl) + \
                  self._get_site_amplification_linear(sites.vs30, C_SR) + \
                  self._get_site_amplification_non_linear(sites.vs30, pga4nl,
                                                          C_SR)
          else:
              mean = self._compute_magnitude_scaling(rup, C) + \
                  self._compute_distance_scaling(rup, dists, C) + \
                  self._get_site_amplification_linear(sites.vs30, C_SR) + \
                  self._get_site_amplification_non_linear(sites.vs30, pga4nl,
                                                          C_SR)
  
          stddevs = self._get_stddevs(C, stddev_types, num_sites=len(sites.vs30))
  
          return mean, stddevs