docs: def get_stddevs(self, C, stddev_shape, stddev_types, sites):
          """
          Returns the standard deviations, with different site standard
          deviation for inferred vs. observed vs30 sites.
          """
          stddevs = []
          tau = C["tau_event"]
          sigma_s = np.zeros(sites.vs30measured.shape, dtype=float)
          sigma_s[sites.vs30measured] += C["sigma_s_obs"]
          sigma_s[np.logical_not(sites.vs30measured)] += C["sigma_s_inf"]
          phi = np.sqrt(C["phi0"] ** 2.0 + sigma_s ** 2.)
          for stddev_type in stddev_types:
              assert stddev_type in self.DEFINED_FOR_STANDARD_DEVIATION_TYPES
              if stddev_type == const.StdDev.TOTAL:
                  stddevs.append(np.sqrt(tau ** 2. + phi ** 2.) +
                                 np.zeros(stddev_shape))
              elif stddev_type == const.StdDev.INTRA_EVENT:
                  stddevs.append(phi + np.zeros(stddev_shape))
              elif stddev_type == const.StdDev.INTER_EVENT:
                  stddevs.append(tau + np.zeros(stddev_shape))
          return stddevs