diff --git a/src/sccala/gen_testdata.py b/src/sccala/gen_testdata.py index 0f3f0a5..7862dde 100644 --- a/src/sccala/gen_testdata.py +++ b/src/sccala/gen_testdata.py @@ -27,6 +27,7 @@ def generate_observed_data( hubble=False, calib=False, h0=70.0, + outl_frac=0.0, ): if rng is None: rng = np.random.default_rng() @@ -36,6 +37,11 @@ def generate_observed_data( col = rng.normal(loc=col_range[0], scale=col_range[1]) ae = np.absolute(rng.normal(loc=ae_range[0], scale=ae_range[1])) + # Optional outlier generation + if outl_frac > 0 and rng.uniform() < outl_frac: + vel = np.absolute(rng.normal(loc=vel_range[0], scale=5000e3)) + col = rng.normal(loc=col_range[0], scale=0.5) + ae = np.absolute(rng.normal(loc=ae_range[0], scale=0.5)) vel_obs = rng.normal(loc=vel, scale=verr) col_obs = rng.normal(loc=col, scale=cerr) ae_obs = rng.normal(loc=ae, scale=aeerr) @@ -120,6 +126,7 @@ def gen_testdata( sig_cut_nom=None, calib_m_cut_nom=None, calib_sig_cut_nom=None, + outl_frac=0.0, ): """ Function generating simulated datasets for standardisation @@ -173,6 +180,8 @@ def gen_testdata( calib_m_cut_nom : float calib_sig_cut_nom : float Parameters for the nominal detection values exported to the data. + outl_frac : float + Fraction of outliers to be generated. Default: 0.0 Returns ------- @@ -214,6 +223,7 @@ def gen_testdata( rng=rng, hubble=hubble, h0=h0, + outl_frac=outl_frac, ) detection_prob = detection_probability(mag, m_cut=m_cut, sigma_cut=sigma_cut) detected = rng.uniform() < detection_prob @@ -302,6 +312,7 @@ def gen_testdata( hubble=hubble, calib=True, h0=h0, + outl_frac=outl_frac, ) detection_prob = detection_probability( mag, m_cut=calib_m_cut, sigma_cut=calib_sigma_cut @@ -463,6 +474,7 @@ def main(args): calib_sigma_cut=args.calib_sigma_cut, calib_m_cut_nom=args.calib_m_cut_nom, calib_sig_cut_nom=args.calib_sig_cut_nom, + outl_frac=args.outl_frac, ) return data @@ -564,6 +576,12 @@ def cli(): type=float, help="Nominal sigma cut for the exported calibrator data. Default: 0.5", ) + parser.add_argument( + "--outl_frac", + type=float, + help="Fraction of outliers to be generated. Default: 0.0", + default=0.0, + ) args = parser.parse_args() diff --git a/src/sccala/sccala_scm.py b/src/sccala/sccala_scm.py index 08ac11e..b8861f9 100644 --- a/src/sccala/sccala_scm.py +++ b/src/sccala/sccala_scm.py @@ -12,13 +12,19 @@ def main(args): blind=args.unblind, blindkey=args.blindkey, ) + if args.outlier and args.classic: + raise ValueError("Outlier model is not available for classic SCM") model = args.model if not args.classic: if model == "hubble": - model = models.HubbleSCM() + model = models.HubbleSCMOutlier() if args.outlier else models.HubbleSCM() elif model == "hubble-free": - model = models.HubbleFreeSCM() + model = ( + models.HubbleFreeSCMOutlier() + if args.outlier + else models.HubbleFreeSCM() + ) elif model == "hubble-nh": model = models.NHHubbleSCM() elif model == "hubble-nh-simple": @@ -52,6 +58,7 @@ def main(args): save_warmup=args.save_warmup, quiet=False, classic=args.classic, + outlier=args.outlier, output_dir=args.output_dir, test_data=args.test_data, selection_effects=args.selection_effects, @@ -144,6 +151,11 @@ def cli(): action="store_true", help="If flag is given, classical SCM is used instead of extended SCM", ) + parser.add_argument( + "--outlier", + action="store_true", + help="If flag is given, outlier model will be used.", + ) parser.add_argument( "--unblind", action="store_false", diff --git a/src/sccala/scmlib/models.py b/src/sccala/scmlib/models.py index dfdaf9a..f7f9053 100644 --- a/src/sccala/scmlib/models.py +++ b/src/sccala/scmlib/models.py @@ -5,6 +5,7 @@ from sccala.utillib.aux import NumpyEncoder +from sccala.utillib.const import VS_INIT, RV_INIT, VTRUE_INIT class SCM_Model: @@ -141,9 +142,9 @@ def __init__(self): def set_initial_conditions(self, init=None): if init is None: self.init = { - "vs": 7500e3, - "rv": 1000e3, - "v_true": [7500e3] * self.data["sn_idx"], + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], } else: self.init = init @@ -156,6 +157,29 @@ def print_results(self, df, blind=True): return +class HubbleFreeSCMOutlier(SCM_Model): + def __init__(self): + self.__name__ = "hubble-free-scm-outlier" + super().__init__() + + def set_initial_conditions(self, init=None): + if init is None: + self.init = { + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], + } + else: + self.init = init + return + + def print_results(self, df, blind=True): + keys = ["alpha", "beta", "gamma", "sigma_int", "Mi", "outl_frac"] + for key in keys: + print("%s = %.2g +/- %.2g" % (key, np.mean(df[key]), np.std(df[key]))) + return + + class HubbleSCM(SCM_Model): def __init__(self): self.__name__ = "hubble-scm" @@ -165,14 +189,14 @@ def __init__(self): def set_initial_conditions(self, init=None): if init is None: self.init = { - "vs": 7500e3, - "rv": 1000e3, - "v_true": [7500e3] * self.data["sn_idx"], - "calib_v_true": [7500e3] * self.data["calib_sn_idx"], + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], + "calib_v_true": [VTRUE_INIT] * self.data["calib_sn_idx"], } for i in range(self.data["num_calib_dset"]): - self.init["calib_vs.%d" % (i + 1)] = 7500e3 - self.init["calib_rv.%d" % (i + 1)] = 1000e3 + self.init["calib_vs.%d" % (i + 1)] = VS_INIT + self.init["calib_rv.%d" % (i + 1)] = RV_INIT else: self.init = init return @@ -187,6 +211,37 @@ def print_results(self, df, blind=True): return +class HubbleSCMOutlier(SCM_Model): + def __init__(self): + self.__name__ = "hubble-scm-outlier" + super().__init__() + self.hubble = True + + def set_initial_conditions(self, init=None): + if init is None: + self.init = { + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], + "calib_v_true": [VTRUE_INIT] * self.data["calib_sn_idx"], + } + for i in range(self.data["num_calib_dset"]): + self.init["calib_vs.%d" % (i + 1)] = VS_INIT + self.init["calib_rv.%d" % (i + 1)] = RV_INIT + else: + self.init = init + return + + def print_results(self, df, blind=True): + if blind: + keys = ["alpha", "beta", "gamma", "sigma_int", "Mi", "outl_frac"] + else: + keys = ["alpha", "beta", "gamma", "sigma_int", "Mi", "H0", "outl_frac"] + for key in keys: + print("%s = %.4g +/- %.4g" % (key, np.mean(df[key]), np.std(df[key]))) + return + + class ClassicNHHubbleFreeSCM(SCM_Model): def __init__(self): self.__name__ = "classic-nh-hubble-free-scm" @@ -240,9 +295,9 @@ def __init__(self): def set_initial_conditions(self, init=None): if init is None: self.init = { - "vs": 7500e3, - "rv": 1000e3, - "v_true": [7500e3] * self.data["sn_idx"], + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], } else: self.init = init @@ -264,14 +319,14 @@ def __init__(self): def set_initial_conditions(self, init=None): if init is None: self.init = { - "vs": 7500e3, - "rv": 1000e3, - "v_true": [7500e3] * self.data["sn_idx"], - "calib_v_true": [7500e3] * self.data["calib_sn_idx"], + "vs": VS_INIT, + "rv": RV_INIT, + "v_true": [VTRUE_INIT] * self.data["sn_idx"], + "calib_v_true": [VTRUE_INIT] * self.data["calib_sn_idx"], } for i in range(self.data["num_calib_dset"]): - self.init["calib_vs.%d" % (i + 1)] = 7500e3 - self.init["calib_rv.%d" % (i + 1)] = 1000e3 + self.init["calib_vs.%d" % (i + 1)] = VS_INIT + self.init["calib_rv.%d" % (i + 1)] = RV_INIT else: self.init = init return diff --git a/src/sccala/scmlib/models/classic-hubble-free-scm.stan b/src/sccala/scmlib/models/classic-hubble-free-scm.stan index 28460d1..68d243c 100644 --- a/src/sccala/scmlib/models/classic-hubble-free-scm.stan +++ b/src/sccala/scmlib/models/classic-hubble-free-scm.stan @@ -25,6 +25,7 @@ parameters { } transformed parameters{ array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; array[sn_idx] real mean; array[sn_idx] real v_mi; real sigma_int; @@ -36,6 +37,14 @@ transformed parameters{ v_mi[i] = (errors[i][1,1] + sigma_int^2) + sigma_cut^2 + (alpha * rv / (vs * log10()))^2 + (beta * rc)^2; } } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } } model { Mi ~ uniform(-30,0); @@ -43,10 +52,10 @@ model { beta ~ uniform(-20,20); log_sigma ~ uniform(-3,0); - vs ~ cauchy(7500e3,1500e3); + vs ~ cauchy(0.75,0.15); cs ~ cauchy(0,0.5); - rv ~ normal(0,1500e3); + rv ~ normal(0,0.15); rc ~ normal(0,0.5); v_true ~ normal(vs,rv); @@ -55,12 +64,9 @@ model { mag_cut ~ normal(m_cut_nom,0.5); sigma_cut ~ normal(sig_cut_nom,0.25); - for (i in 1:sn_idx) { - target += multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i]]', errors[i] + [[sigma_int^2, 0, 0], [0, 0, 0], [0, 0, 0]]); - if (use_selection != 0) { - target += normal_lcdf(mag_cut | obs[i][1], sigma_cut) - - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + outl_frac ~ lognormal(-3,0.25); - } + for (i in 1:sn_idx) { + target += sn_log_like[i]; } } diff --git a/src/sccala/scmlib/models/classic-hubble-scm.stan b/src/sccala/scmlib/models/classic-hubble-scm.stan index fca408b..f942b1d 100644 --- a/src/sccala/scmlib/models/classic-hubble-scm.stan +++ b/src/sccala/scmlib/models/classic-hubble-scm.stan @@ -28,10 +28,10 @@ parameters { real cs; // Mean of latent color real rv; // Dispersion of latent velocity real rc; // Dispersion of latent color - array[num_calib_dset] real calib_vs; // Mean of latent velocity - array[num_calib_dset] real calib_cs; // Mean of latent color - array[num_calib_dset] real calib_rv; // Dispersion of latent velocity - array[num_calib_dset] real calib_rc; // Dispersion of latent color + //array[num_calib_dset] real calib_vs; // Mean of latent velocity + //array[num_calib_dset] real calib_cs; // Mean of latent color + //array[num_calib_dset] real calib_rv; // Dispersion of latent velocity + //array[num_calib_dset] real calib_rc; // Dispersion of latent color array[sn_idx] real v_true; // Modeled latent velocities (cannot be negative) array[sn_idx] real c_true; // Modeled latent color array[calib_sn_idx] real calib_v_true; // Modeled latent velocities (cannot be negative) @@ -43,9 +43,11 @@ parameters { } transformed parameters{ array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; array[sn_idx] real mean; array[sn_idx] real v_mi; array[calib_sn_idx] real calib_mag_true; + array[calib_sn_idx] real calib_sn_log_like; // array[calib_sn_idx] real calib_mean; // array[calib_sn_idx] real calib_v_mi; real sigma_int; @@ -68,6 +70,21 @@ transformed parameters{ // calib_v_mi[i] = (calib_errors[i][1] + calib_sigma_int[calib_dset_idx[i]]^2) + calib_sigma_cut[calib_dset_idx[i]]^2 + (alpha * calib_rv[calib_dset_idx[i]] / (calib_vs[calib_dset_idx[i]] * log10()))^2 + (beta * calib_rc[calib_dset_idx[i]])^2; // } } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } + for (i in 1:calib_sn_idx) { + calib_sn_log_like[i] = normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i]]', sqrt(calib_errors[i] + [calib_sigma_int[calib_dset_idx[i]]^2, 0, 0]')); + // if (use_selection != 0) { + // calib_sn_log_like[i] += normal_lcdf(calib_mag_cut[calib_dset_idx[i]] | calib_obs[i][1], calib_sigma_cut[calib_dset_idx[i]]) + // - log(normal_cdf(calib_mag_cut[calib_dset_idx[i]] | calib_mean[i], sqrt(calib_v_mi[i])) + 0.0001); + // } + } } model { H0 ~ uniform(0,200); @@ -79,26 +96,28 @@ model { calib_log_sigma[i] ~ uniform(-3,0); } - vs ~ cauchy(7500e3,1500e3); + vs ~ cauchy(0.75,0.15); cs ~ cauchy(0,0.5); - rv ~ normal(0,1500e3); + rv ~ normal(0,0.15); rc ~ normal(0,0.5); v_true ~ normal(vs,rv); c_true ~ normal(cs,rc); - for (i in 1:num_calib_dset) { - calib_vs[i] ~ cauchy(7500e3,1500e3); - calib_cs[i] ~ cauchy(0,0.5); + //for (i in 1:num_calib_dset) { + // calib_vs[i] ~ cauchy(0.75,0.15); + // calib_cs[i] ~ cauchy(0,0.5); - calib_rv[i] ~ normal(0,1500e3); - calib_rc[i] ~ normal(0,0.5); - } + // calib_rv[i] ~ normal(0,0.15); + // calib_rc[i] ~ normal(0,0.5); + //} for (i in 1:calib_sn_idx) { - calib_v_true[i] ~ normal(calib_vs[calib_dset_idx[i]],calib_rv[calib_dset_idx[i]]); - calib_c_true[i] ~ normal(calib_cs[calib_dset_idx[i]],calib_rc[calib_dset_idx[i]]); + //calib_v_true[i] ~ normal(calib_vs[calib_dset_idx[i]],calib_rv[calib_dset_idx[i]]); + //calib_c_true[i] ~ normal(calib_cs[calib_dset_idx[i]],calib_rc[calib_dset_idx[i]]); + calib_v_true[i] ~ normal(vs,rv); + calib_c_true[i] ~ normal(cs,rc); } mag_cut ~ normal(m_cut_nom,0.5); @@ -110,17 +129,9 @@ model { // } for (i in 1:sn_idx) { - target += multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i]]', errors[i] + [[sigma_int^2, 0, 0], [0, 0, 0], [0, 0, 0]]); - if (use_selection != 0) { - target += normal_lcdf(mag_cut | obs[i][1], sigma_cut) - - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); - } + target += sn_log_like[i]; } for (i in 1:calib_sn_idx) { - target += normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i]]', sqrt(calib_errors[i] + [calib_sigma_int[calib_dset_idx[i]]^2, 0, 0]')); - // if (use_selection != 0) { - // target += normal_lcdf(calib_mag_cut[calib_dset_idx[i]] | calib_obs[i][1], calib_sigma_cut[calib_dset_idx[i]]) - // - log(normal_cdf(calib_mag_cut[calib_dset_idx[i]] | calib_mean[i], sqrt(calib_v_mi[i])) + 0.0001); - // } + target += calib_sn_log_like[i]; } } diff --git a/src/sccala/scmlib/models/hubble-free-scm-outlier.stan b/src/sccala/scmlib/models/hubble-free-scm-outlier.stan new file mode 100644 index 0000000..6404254 --- /dev/null +++ b/src/sccala/scmlib/models/hubble-free-scm-outlier.stan @@ -0,0 +1,94 @@ +data { + int sn_idx; + array[sn_idx] vector[4] obs; // Observed SN properties + array[sn_idx] matrix[4,4] errors; // Associated uncertaintes (measurement, statistical, systematic) + real vel_avg; // Normalisation constants + real col_avg; + real ae_avg; + array[sn_idx] real log_dist_mod; // Pre-computed, redshift dependent, Hubble-free distance moduli + real m_cut_nom; // Nominal magnitude cut for selection effect calculation + real sig_cut_nom; // Nominal uncertainty of the magnitude cut + int use_selection; // Flag to use selection effect. If 0, the selection effect is not used + vector[4] outl_width; // Width of the outlier distribution +} +parameters { + real Mi; // Absolute Hubble-free Magnitude + real alpha; // Velocity correction strength + real beta; // Color correction strength + real gamma; // a/e correction strength + real log_sigma; // Unexplained intrinsic scatter + real vs; // Mean of latent velocity + real cs; // Mean of latent color + real as; // Mean of latent a/e + real rv; // Dispersion of latent velocity + real rc; // Dispersion of latent color + real ra; // Dispersion of latent a/e + array[sn_idx] real v_true; // Modeled latent velocities (cannot be negative) + array[sn_idx] real c_true; // Modeled latent color + array[sn_idx] real a_true; // Modeled latent a/e (cannot be negative) + real mag_cut; // Magnitude cut for selection effect calculation + real sigma_cut; // Uncertainty of the magnitude cut + real log_outl_frac; +} +transformed parameters{ + array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; + array[sn_idx] real mean; + array[sn_idx] real v_mi; + real sigma_int; + sigma_int = 10 ^ log_sigma; + real outl_frac; + for (i in 1:sn_idx) { + mag_true[i] = Mi - alpha * log10(v_true[i] / vel_avg) + beta * (c_true[i] - col_avg) + gamma * (a_true[i] - ae_avg) + 5 * log_dist_mod[i]; + if (use_selection != 0) { + mean[i] = Mi - alpha * log10(vs / vel_avg) + beta * (cs - col_avg) + gamma * (as - ae_avg) + 5 * log_dist_mod[i]; + v_mi[i] = (errors[i][1,1] + sigma_int^2) + sigma_cut^2 + (alpha * rv / (vs * log10()))^2 + (beta * rc)^2 + (gamma * ra)^2; + } + } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } + outl_frac = 10 ^ log_outl_frac; +} +model { + Mi ~ uniform(-30,0); + alpha ~ uniform(-20,20); + beta ~ uniform(-20,20); + gamma ~ uniform(-20,20); + log_sigma ~ uniform(-3,0); + + vs ~ cauchy(0.75,0.15); + cs ~ cauchy(0,0.5); + as ~ cauchy(0.5,0.5); + + rv ~ normal(0,0.15); + rc ~ normal(0,0.5); + ra ~ normal(0,0.5); + + v_true ~ normal(vs,rv); + c_true ~ normal(cs,rc); + a_true ~ normal(as,ra); + + mag_cut ~ normal(m_cut_nom,0.5); + sigma_cut ~ normal(sig_cut_nom,0.25); + + log_outl_frac ~ uniform(-5,-0.2); + + for (i in 1:sn_idx) { + target += log_sum_exp( + (log(1 - outl_frac) + sn_log_like[i]), + (log(outl_frac) + multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', diag_matrix(outl_width))) + ); + } +} +generated quantities { + array[sn_idx] real outl_log_like; + for (i in 1:sn_idx) { + outl_log_like[i] = log(1 - outl_frac) + sn_log_like[i] - (log(outl_frac) + multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', diag_matrix(outl_width))); + } +} diff --git a/src/sccala/scmlib/models/hubble-free-scm.stan b/src/sccala/scmlib/models/hubble-free-scm.stan index 6f8c6bc..701a4d9 100644 --- a/src/sccala/scmlib/models/hubble-free-scm.stan +++ b/src/sccala/scmlib/models/hubble-free-scm.stan @@ -30,6 +30,7 @@ parameters { } transformed parameters{ array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; array[sn_idx] real mean; array[sn_idx] real v_mi; real sigma_int; @@ -41,6 +42,14 @@ transformed parameters{ v_mi[i] = (errors[i][1,1] + sigma_int^2) + sigma_cut^2 + (alpha * rv / (vs * log10()))^2 + (beta * rc)^2 + (gamma * ra)^2; } } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } } model { Mi ~ uniform(-30,0); @@ -49,11 +58,11 @@ model { gamma ~ uniform(-20,20); log_sigma ~ uniform(-3,0); - vs ~ cauchy(7500e3,1500e3); + vs ~ cauchy(0.75,0.15); cs ~ cauchy(0,0.5); as ~ cauchy(0.5,0.5); - rv ~ normal(0,1500e3); + rv ~ normal(0,0.15); rc ~ normal(0,0.5); ra ~ normal(0,0.5); @@ -65,10 +74,6 @@ model { sigma_cut ~ normal(sig_cut_nom,0.25); for (i in 1:sn_idx) { - target += multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + [[sigma_int^2, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]); - if (use_selection != 0) { - target += normal_lcdf(mag_cut | obs[i][1], sigma_cut) - - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); - } + target += sn_log_like[i]; } } diff --git a/src/sccala/scmlib/models/hubble-scm-outlier.stan b/src/sccala/scmlib/models/hubble-scm-outlier.stan new file mode 100644 index 0000000..4be8252 --- /dev/null +++ b/src/sccala/scmlib/models/hubble-scm-outlier.stan @@ -0,0 +1,175 @@ +data { + int sn_idx; + array[sn_idx] vector[4] obs; // Observed SN properties + array[sn_idx] matrix[4,4] errors; // Associated uncertaintes (measurement, statistical, systematic) + real vel_avg; // Normalisation constans + real col_avg; + real ae_avg; + array[sn_idx] real log_dist_mod; // Pre-computed, redshift dependent, Hubble-free distance moduli + int calib_sn_idx; + array[calib_sn_idx] vector[4] calib_obs; // Observed SN properties + array[calib_sn_idx] vector[4] calib_errors; // Associated uncertaintes (measurement, statistical, systematic) + array[calib_sn_idx] real calib_dist_mod; // Distance moduli of calibrators + array[calib_sn_idx] int calib_dset_idx; // Index of the calibrator dataset + int num_calib_dset; // Number of calibrator datasets + real m_cut_nom; // Nominal magnitude cut for selection effect calculation + real sig_cut_nom; // Nominal uncertainty of the magnitude cut + // array[num_calib_dset] real calib_m_cut_nom; // Nominal magnitude cut for selection effect calculation + // array[num_calib_dset] real calib_sig_cut_nom; // Nominal uncertainty of the magnitude cut + int use_selection; // Flag to use selection effect. If 0, the selection effect is not used + vector[4] outl_width; // Width of the outlier distribution +} +parameters { + real H0; // Hubble constant + real Mi; // Absolute Hubble-free Magnitude + real alpha; // Velocity correction strength + real beta; // Color correction strength + real gamma; // a/e correction strength + real log_sigma; // Unexplained intrinsic scatter + array[num_calib_dset] real calib_log_sigma; // Unexplained intrinsic scatter + real vs; // Mean of latent velocity + real cs; // Mean of latent color + real as; // Mean of latent a/e + real rv; // Dispersion of latent velocity + real rc; // Dispersion of latent color + real ra; // Dispersion of latent a/e + //array[num_calib_dset] real calib_vs; // Mean of latent velocity + //array[num_calib_dset] real calib_cs; // Mean of latent color + //array[num_calib_dset] real calib_as; // Mean of latent a/e + //array[num_calib_dset] real calib_rv; // Dispersion of latent velocity + //array[num_calib_dset] real calib_rc; // Dispersion of latent color + //array[num_calib_dset] real calib_ra; // Dispersion of latent a/e + array[sn_idx] real v_true; // Modeled latent velocities (cannot be negative) + array[sn_idx] real c_true; // Modeled latent color + array[sn_idx] real a_true; // Modeled latent a/e (cannot be negative) + array[calib_sn_idx] real calib_v_true; // Modeled latent velocities (cannot be negative) + array[calib_sn_idx] real calib_c_true; // Modeled latent color + array[calib_sn_idx] real calib_a_true; // Modeled latent a/e (cannot be negative) + real mag_cut; // Magnitude cut for selection effect calculation + real sigma_cut; // Uncertainty of the magnitude cut + // array[num_calib_dset] real calib_mag_cut; // Magnitude cut for selection effect calculation + // array[num_calib_dset] real calib_sigma_cut; // Uncertainty of the magnitude cut + real log_outl_frac; +} +transformed parameters{ + array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; + array[sn_idx] real mean; + array[sn_idx] real v_mi; + array[calib_sn_idx] real calib_mag_true; + array[calib_sn_idx] real calib_sn_log_like; + // array[calib_sn_idx] real calib_mean; + // array[calib_sn_idx] real calib_v_mi; + real sigma_int; + array[num_calib_dset] real calib_sigma_int; + real outl_frac; + sigma_int = 10 ^ log_sigma; + for (i in 1:num_calib_dset) { + calib_sigma_int[i] = 10 ^ calib_log_sigma[i]; + } + for (i in 1:sn_idx) { + mag_true[i] = Mi - 5 * log10(H0) + 25 - alpha * log10(v_true[i] / vel_avg) + beta * (c_true[i] - col_avg) + gamma * (a_true[i] - ae_avg) + 5 * log_dist_mod[i]; + if (use_selection != 0) { + mean[i] = Mi - 5 * log10(H0) + 25 - alpha * log10(vs / vel_avg) + beta * (cs - col_avg) + gamma * (as - ae_avg) + 5 * log_dist_mod[i]; + v_mi[i] = (errors[i][1,1] + sigma_int^2) + sigma_cut^2 + (alpha * rv / (vs * log10()))^2 + (beta * rc)^2 + (gamma * ra)^2; + } + } + for (i in 1:calib_sn_idx) { + calib_mag_true[i] = Mi - alpha * log10(calib_v_true[i] / vel_avg) + beta * (calib_c_true[i] - col_avg) + gamma * (calib_a_true[i] - ae_avg) + calib_dist_mod[i]; + // if (use_selection != 0) { + // calib_mean[i] = Mi - alpha * log10(calib_vs[calib_dset_idx[i]] / vel_avg) + beta * (calib_cs[calib_dset_idx[i]] - col_avg) + gamma * (calib_as[calib_dset_idx[i]] - ae_avg) + calib_dist_mod[i]; + // calib_v_mi[i] = (calib_errors[i][1] + calib_sigma_int[calib_dset_idx[i]]^2) + calib_sigma_cut[calib_dset_idx[i]]^2 + (alpha * calib_rv[calib_dset_idx[i]] / (calib_vs[calib_dset_idx[i]] * log10()))^2 + (beta * calib_rc[calib_dset_idx[i]])^2 + (gamma * calib_ra[calib_dset_idx[i]])^2; + // } + } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } + for (i in 1:calib_sn_idx) { + calib_sn_log_like[i] = normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i], calib_a_true[i]]', sqrt(calib_errors[i] + [calib_sigma_int[calib_dset_idx[i]]^2, 0, 0, 0]')); + // if (use_selection != 0) { + // calib_sn_log_like[i] += normal_lcdf(calib_mag_cut[calib_dset_idx[i]] | calib_obs[i][1], calib_sigma_cut[calib_dset_idx[i]]) + // - log(normal_cdf(calib_mag_cut[calib_dset_idx[i]] | calib_mean[i], sqrt(calib_v_mi[i])) + 0.0001); + // } + } + outl_frac = 10 ^ log_outl_frac; +} +model { + H0 ~ uniform(0,200); + Mi ~ uniform(-30,0); + alpha ~ uniform(-20,20); + beta ~ uniform(-20,20); + gamma ~ uniform(-20,20); + log_sigma ~ uniform(-3,0); + for (i in 1:num_calib_dset) { + calib_log_sigma[i] ~ uniform(-3,0); + } + + vs ~ cauchy(0.75,0.15); + cs ~ cauchy(0,0.5); + as ~ cauchy(0.5,0.5); + + rv ~ normal(0,0.15); + rc ~ normal(0,0.5); + ra ~ normal(0,0.5); + + v_true ~ normal(vs,rv); + c_true ~ normal(cs,rc); + a_true ~ normal(as,ra); + + //for (i in 1:num_calib_dset) { + // calib_vs[i] ~ cauchy(0.75,0.15); + // calib_cs[i] ~ cauchy(0,0.5); + // calib_as[i] ~ cauchy(0.5,0.5); + + // calib_rv[i] ~ normal(0,1.5); + // calib_rc[i] ~ normal(0,0.5); + // calib_ra[i] ~ normal(0,0.5); + //} + + for (i in 1:calib_sn_idx) { + //calib_v_true[i] ~ normal(calib_vs[calib_dset_idx[i]],calib_rv[calib_dset_idx[i]]); + //calib_c_true[i] ~ normal(calib_cs[calib_dset_idx[i]],calib_rc[calib_dset_idx[i]]); + //calib_a_true[i] ~ normal(calib_as[calib_dset_idx[i]],calib_ra[calib_dset_idx[i]]); + calib_v_true[i] ~ normal(vs,rv); + calib_c_true[i] ~ normal(cs,rc); + calib_a_true[i] ~ normal(as,ra); + } + + mag_cut ~ normal(m_cut_nom,0.5); + sigma_cut ~ normal(sig_cut_nom,0.25); + + log_outl_frac ~ uniform(-5,-0.2); + + // for (i in 1:calib_sn_idx) { + // calib_mag_cut[calib_dset_idx[i]] ~ normal(calib_m_cut_nom[calib_dset_idx[i]],0.5); + // calib_sigma_cut[calib_dset_idx[i]] ~ normal(calib_sig_cut_nom[calib_dset_idx[i]],0.25); + // } + + for (i in 1:sn_idx) { + target += log_sum_exp( + (log(1 - outl_frac) + sn_log_like[i]), + (log(outl_frac) + multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', diag_matrix(outl_width))) + ); + } + for (i in 1:calib_sn_idx) { + target += log_sum_exp( + (log(1 - outl_frac) + calib_sn_log_like[i]), + (log(outl_frac) + multi_normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i], calib_a_true[i]]', diag_matrix(outl_width))) + ); + } +} +generated quantities { + array[sn_idx] real outl_log_like; + array[calib_sn_idx] real calib_outl_log_like; + for (i in 1:sn_idx) { + outl_log_like[i] = log(1 - outl_frac) + sn_log_like[i] - (log(outl_frac) + multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', diag_matrix(outl_width))); + } + for (i in 1:calib_sn_idx) { + calib_outl_log_like[i] = log(1 - outl_frac) + calib_sn_log_like[i] - (log(outl_frac) + multi_normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i], calib_a_true[i]]', diag_matrix(outl_width))); + } +} diff --git a/src/sccala/scmlib/models/hubble-scm.stan b/src/sccala/scmlib/models/hubble-scm.stan index 18b46e8..a6a8790 100644 --- a/src/sccala/scmlib/models/hubble-scm.stan +++ b/src/sccala/scmlib/models/hubble-scm.stan @@ -32,12 +32,12 @@ parameters { real rv; // Dispersion of latent velocity real rc; // Dispersion of latent color real ra; // Dispersion of latent a/e - array[num_calib_dset] real calib_vs; // Mean of latent velocity - array[num_calib_dset] real calib_cs; // Mean of latent color - array[num_calib_dset] real calib_as; // Mean of latent a/e - array[num_calib_dset] real calib_rv; // Dispersion of latent velocity - array[num_calib_dset] real calib_rc; // Dispersion of latent color - array[num_calib_dset] real calib_ra; // Dispersion of latent a/e + //array[num_calib_dset] real calib_vs; // Mean of latent velocity + //array[num_calib_dset] real calib_cs; // Mean of latent color + //array[num_calib_dset] real calib_as; // Mean of latent a/e + //array[num_calib_dset] real calib_rv; // Dispersion of latent velocity + //array[num_calib_dset] real calib_rc; // Dispersion of latent color + //array[num_calib_dset] real calib_ra; // Dispersion of latent a/e array[sn_idx] real v_true; // Modeled latent velocities (cannot be negative) array[sn_idx] real c_true; // Modeled latent color array[sn_idx] real a_true; // Modeled latent a/e (cannot be negative) @@ -51,9 +51,11 @@ parameters { } transformed parameters{ array[sn_idx] real mag_true; + array[sn_idx] real sn_log_like; array[sn_idx] real mean; array[sn_idx] real v_mi; array[calib_sn_idx] real calib_mag_true; + array[calib_sn_idx] real calib_sn_log_like; // array[calib_sn_idx] real calib_mean; // array[calib_sn_idx] real calib_v_mi; real sigma_int; @@ -76,6 +78,21 @@ transformed parameters{ // calib_v_mi[i] = (calib_errors[i][1] + calib_sigma_int[calib_dset_idx[i]]^2) + calib_sigma_cut[calib_dset_idx[i]]^2 + (alpha * calib_rv[calib_dset_idx[i]] / (calib_vs[calib_dset_idx[i]] * log10()))^2 + (beta * calib_rc[calib_dset_idx[i]])^2 + (gamma * calib_ra[calib_dset_idx[i]])^2; // } } + + for (i in 1:sn_idx) { + sn_log_like[i] = multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + diag_matrix([sigma_int^2, 0, 0, 0]')); + if (use_selection != 0) { + sn_log_like[i] += normal_lcdf(mag_cut | obs[i][1], sigma_cut) + - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); + } + } + for (i in 1:calib_sn_idx) { + calib_sn_log_like[i] = normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i], calib_a_true[i]]', sqrt(calib_errors[i] + [calib_sigma_int[calib_dset_idx[i]]^2, 0, 0, 0]')); + // if (use_selection != 0) { + // calib_sn_log_like[i] += normal_lcdf(calib_mag_cut[calib_dset_idx[i]] | calib_obs[i][1], calib_sigma_cut[calib_dset_idx[i]]) + // - log(normal_cdf(calib_mag_cut[calib_dset_idx[i]] | calib_mean[i], sqrt(calib_v_mi[i])) + 0.0001); + // } + } } model { H0 ~ uniform(0,200); @@ -88,11 +105,11 @@ model { calib_log_sigma[i] ~ uniform(-3,0); } - vs ~ cauchy(7500e3,1500e3); + vs ~ cauchy(0.75,0.15); cs ~ cauchy(0,0.5); as ~ cauchy(0.5,0.5); - rv ~ normal(0,1500e3); + rv ~ normal(0,0.15); rc ~ normal(0,0.5); ra ~ normal(0,0.5); @@ -100,20 +117,23 @@ model { c_true ~ normal(cs,rc); a_true ~ normal(as,ra); - for (i in 1:num_calib_dset) { - calib_vs[i] ~ cauchy(7500e3,1500e3); - calib_cs[i] ~ cauchy(0,0.5); - calib_as[i] ~ cauchy(0.5,0.5); + //for (i in 1:num_calib_dset) { + // calib_vs[i] ~ cauchy(0.75,0.15); + // calib_cs[i] ~ cauchy(0,0.5); + // calib_as[i] ~ cauchy(0.5,0.5); - calib_rv[i] ~ normal(0,1500e3); - calib_rc[i] ~ normal(0,0.5); - calib_ra[i] ~ normal(0,0.5); - } + // calib_rv[i] ~ normal(0,0.15); + // calib_rc[i] ~ normal(0,0.5); + // calib_ra[i] ~ normal(0,0.5); + //} for (i in 1:calib_sn_idx) { - calib_v_true[i] ~ normal(calib_vs[calib_dset_idx[i]],calib_rv[calib_dset_idx[i]]); - calib_c_true[i] ~ normal(calib_cs[calib_dset_idx[i]],calib_rc[calib_dset_idx[i]]); - calib_a_true[i] ~ normal(calib_as[calib_dset_idx[i]],calib_ra[calib_dset_idx[i]]); + //calib_v_true[i] ~ normal(calib_vs[calib_dset_idx[i]],calib_rv[calib_dset_idx[i]]); + //calib_c_true[i] ~ normal(calib_cs[calib_dset_idx[i]],calib_rc[calib_dset_idx[i]]); + //calib_a_true[i] ~ normal(calib_as[calib_dset_idx[i]],calib_ra[calib_dset_idx[i]]); + calib_v_true[i] ~ normal(vs,rv); + calib_c_true[i] ~ normal(cs,rc); + calib_a_true[i] ~ normal(as,ra); } mag_cut ~ normal(m_cut_nom,0.5); @@ -125,17 +145,9 @@ model { // } for (i in 1:sn_idx) { - target += multi_normal_lpdf(obs[i] | [mag_true[i], v_true[i], c_true[i], a_true[i]]', errors[i] + [[sigma_int^2, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0], [0, 0, 0, 0]]); - if (use_selection != 0) { - target += normal_lcdf(mag_cut | obs[i][1], sigma_cut) - - log(normal_cdf(mag_cut | mean[i], sqrt(v_mi[i])) + 0.0001); - } + target += sn_log_like[i]; } for (i in 1:calib_sn_idx) { - target += normal_lpdf(calib_obs[i] | [calib_mag_true[i], calib_v_true[i], calib_c_true[i], calib_a_true[i]]', sqrt(calib_errors[i] + [calib_sigma_int[calib_dset_idx[i]]^2, 0, 0, 0]')); - // if (use_selection != 0) { - // target += normal_lcdf(calib_mag_cut[calib_dset_idx[i]] | calib_obs[i][1], calib_sigma_cut[calib_dset_idx[i]]) - // - log(normal_cdf(calib_mag_cut[calib_dset_idx[i]] | calib_mean[i], sqrt(calib_v_mi[i])) + 0.0001); - // } + target += calib_sn_log_like[i]; } } diff --git a/src/sccala/scmlib/sccala.py b/src/sccala/scmlib/sccala.py index 546bbf0..7450bb5 100644 --- a/src/sccala/scmlib/sccala.py +++ b/src/sccala/scmlib/sccala.py @@ -12,7 +12,7 @@ from sccala.scmlib.models import SCM_Model from sccala.utillib.aux import distmod_kin, quantile, split_list, nullify_output -from sccala.utillib.const import C_LIGHT, PV_UNCERTAINTY +from sccala.utillib.const import C_LIGHT, PV_UNCERTAINTY, VEL_NORM class SccalaSCM: @@ -61,8 +61,8 @@ def __init__(self, file, calib="CALIB", blind=True, blindkey="HUBBLE"): self.col = df[df["dataset"].isin(datasets)]["col"].to_numpy() self.col_err = df[df["dataset"].isin(datasets)]["col_err"].to_numpy() - self.vel = df[df["dataset"].isin(datasets)]["vel"].to_numpy() - self.vel_err = df[df["dataset"].isin(datasets)]["vel_err"].to_numpy() + self.vel = df[df["dataset"].isin(datasets)]["vel"].to_numpy() / VEL_NORM + self.vel_err = df[df["dataset"].isin(datasets)]["vel_err"].to_numpy() / VEL_NORM self.ae = df[df["dataset"].isin(datasets)]["ae"].to_numpy() self.ae_err = df[df["dataset"].isin(datasets)]["ae_err"].to_numpy() @@ -71,7 +71,7 @@ def __init__(self, file, calib="CALIB", blind=True, blindkey="HUBBLE"): self.red_err = df[df["dataset"].isin(datasets)]["red_err"].to_numpy() self.mag_sys = df[df["dataset"].isin(datasets)]["mag_sys"].to_numpy() - self.v_sys = df[df["dataset"].isin(datasets)]["vel_sys"].to_numpy() + self.v_sys = df[df["dataset"].isin(datasets)]["vel_sys"].to_numpy() / VEL_NORM self.c_sys = df[df["dataset"].isin(datasets)]["col_sys"].to_numpy() self.ae_sys = df[df["dataset"].isin(datasets)]["ae_sys"].to_numpy() @@ -93,10 +93,12 @@ def __init__(self, file, calib="CALIB", blind=True, blindkey="HUBBLE"): "col_err" ].to_numpy() - self.calib_vel = df[df["dataset"].isin(calib_datasets)]["vel"].to_numpy() - self.calib_vel_err = df[df["dataset"].isin(calib_datasets)][ - "vel_err" - ].to_numpy() + self.calib_vel = ( + df[df["dataset"].isin(calib_datasets)]["vel"].to_numpy() / VEL_NORM + ) + self.calib_vel_err = ( + df[df["dataset"].isin(calib_datasets)]["vel_err"].to_numpy() / VEL_NORM + ) self.calib_ae = df[df["dataset"].isin(calib_datasets)]["ae"].to_numpy() self.calib_ae_err = df[df["dataset"].isin(calib_datasets)][ @@ -111,9 +113,9 @@ def __init__(self, file, calib="CALIB", blind=True, blindkey="HUBBLE"): self.calib_mag_sys = df[df["dataset"].isin(calib_datasets)][ "mag_sys" ].to_numpy() - self.calib_v_sys = df[df["dataset"].isin(calib_datasets)][ - "vel_sys" - ].to_numpy() + self.calib_v_sys = ( + df[df["dataset"].isin(calib_datasets)]["vel_sys"].to_numpy() / VEL_NORM + ) self.calib_c_sys = df[df["dataset"].isin(calib_datasets)][ "col_sys" ].to_numpy() @@ -263,6 +265,8 @@ def sample( quiet=False, init=None, classic=False, + outlier=False, + outl_width=[1.0, 1.0, 1.0, 1.0], output_dir=None, test_data=False, selection_effects=True, @@ -297,6 +301,11 @@ def sample( classic : bool Switches classic mode on if True. In classic mode, a/e input is ignored. + outlier : bool + If True, outlier detection is enabled. Incompatible with + classic mode. Default: False + outl_width : list + Width of the outlier detection. Default: [1.0, 1.0, 1.0, 1.0] output_dir : str Directory where temporary STAN files will be stored. Default: None test_data : bool @@ -348,7 +357,7 @@ def sample( model.data["use_selection"] = 0 if test_data: - model.data["vel_avg"] = 7100e3 + model.data["vel_avg"] = 7100e3 / VEL_NORM model.data["col_avg"] = 0.5 if not classic: model.data["ae_avg"] = 0.31 @@ -414,6 +423,14 @@ def sample( # model.data["calib_m_cut_nom"] = np.zeros(n_calib_dset) # model.data["calib_sig_cut_nom"] = np.zeros(n_calib_dset) + model.data["vel_avg"] = np.mean(np.concatenate([self.vel, self.calib_vel])) + model.data["col_avg"] = np.mean(np.concatenate([self.col, self.calib_col])) + if not classic: + model.data["ae_avg"] = np.mean(np.concatenate([self.ae, self.calib_ae])) + + if outlier: + model.data["outl_width"] = np.array(outl_width) + model.set_initial_conditions(init) if log_dir is not None: @@ -1048,6 +1065,10 @@ def cornerplot(self, save=None, classic=False): paramnames.append(r"$\sigma_\mathrm{cut}$") keys.append("sigma_cut") ndim += 1 + if "outl_frac" in all_keys: + paramnames.append(r"$f_\mathrm{outl}$") + keys.append("outl_frac") + ndim += 1 posterior = self.posterior[keys].to_numpy() diff --git a/src/sccala/utillib/const.py b/src/sccala/utillib/const.py index 7f9a896..aadc0a9 100644 --- a/src/sccala/utillib/const.py +++ b/src/sccala/utillib/const.py @@ -2,3 +2,7 @@ H_ERG = 6.62607015e-27 # Plancks constant (erg.s) C_AA = 299792458 * 1.0e10 # in AA/s PV_UNCERTAINTY = 150.0 # Peculiar velocity uncertainty in km/s +VEL_NORM = 10000.0e3 # Velocity normalization in m/s +VS_INIT = 0.75 # Initial guess for the velocity scale in 10000 km/s +RV_INIT = 0.1 # Initial guess for the RV in 10000 km/s +VTRUE_INIT = 0.75 # Initial guess for the true velocity in 10000 km/s