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.
  
          Implement equation 1, page 20.
          """
          # compute median PGA on rock, needed to compute non-linear site
          # amplification
          C_pga = self.COEFFS[PGA()]
          median_pga = np.exp(
              self._compute_mean(C_pga, rup.mag, dists, rup.rake)
          )
  
          # compute full mean value by adding nonlinear site amplification terms
          C = self.COEFFS[imt]
          mean = (self._compute_mean(C, rup.mag, dists, rup.rake) +
                  self._compute_non_linear_term(C, median_pga, sites))
  
          stddevs = self._get_stddevs(C, stddev_types, num_sites=sites.vs30.size)
  
          return mean + self.adjustment_factor, stddevs