code: def get_mean_and_stddevs(self, sites, rup, dists, imt, stddev_types):
          # pylint: disable=too-many-arguments
          
  
          # obtain coefficients for required intensity measure type
          coeffs = self.COEFFS_BEDROCK[imt].copy()
  
          # obtain site-class specific coefficients
          a_1, a_2, sigma_site = self._get_site_coeffs(sites, imt)
          coeffs.update({'a1': a_1, 'a2': a_2, 'sigma_site': sigma_site})
  
          # compute bedrock motion, equation (8)
          ln_mean = (self._compute_magnitude_terms(rup, coeffs) +
                     self._compute_distance_terms(dists, coeffs))
  
          # adjust for site class, equation (10)
          ln_mean += self._compute_site_amplification(ln_mean, coeffs)
          # No need to convert to g since "In [equation (8)], y_br = (SA/g)"
  
          ln_stddevs = self._get_stddevs(coeffs, stddev_types)
  
          return ln_mean, [ln_stddevs]
docs: See :meth:`superclass method
          <.base.GroundShakingIntensityModel.get_mean_and_stddevs>`
          for specification of input and result values.
  
          Implements the following equations:
  
          Equation (8) on p. 203 for the bedrock ground motion:
  
          ``ln(y_br) = c1 + c2*(M - 6) + c3*(M - 6)**2 - lnR - c4*R + ln(ε_br)``
  
          Equation (9) on p. 207 gives the site amplification factor:
  
          ``ln(F_s) = a1*y_br + a2 + ln(δ_site)``
  
          Equation (10) on p. 207 for the ground motion at a given site:
  
          ``y_site = y_br*F_s``
  
          Equation (11) on p. 207 for total standard error at a given site:
  
          ``σ{ln(ε_site)} = sqrt(σ{ln(ε_br)}**2 + σ{ln(δ_site)}**2)``