"""Step 16: robustness of the raw (pre-vetting) outlier statistics to analysis choices, and an empirical look at the Cosmicflows-4 curation. For each variant we refit the relation and report the sample size, the robust scatter, and the number of galaxies below each threshold. (Vetting is not repeated; the point is to show that the size of the faint tail is not created by a particular cut or correction.) The running-median variants replace the clipped polynomial by a 15-bin running median (fit-form check; the full vetted limit for that variant is in limit_variants.csv). CF4 curation: CF4 took its ALFALFA galaxies from the same HI catalogue as our "gold" sample, then applied its own cuts (type, inclination > 45 deg by eye, "no gross anomaly", photometry quality). The fraction of ordinary (catalogue-clean, |D| < 0.3 dex) ALFALFA gold galaxies that ended up in our CF4 sample, as a function of apparent magnitude and surface brightness, shows how strongly that curation disfavours faint / diffuse discs. It is a measurement on NATURAL faint galaxies, not on dimmed ones. Outputs: data/robustness.csv, data/cf4_membership.csv """ import os import numpy as np import pandas as pd from common import DATA, C_KMS, robust_polyfit, mad_sigma, running_median_fit, clopper_pearson X0 = 2.2 def stats(name, logV, bands, deg=2, extra=None): res = [] for y in bands: m = np.isfinite(logV) & np.isfinite(y) if deg == "rm": # 15-bin running median res.append(y - running_median_fit(logV[m], y[m])(logV)) continue p, _ = robust_polyfit(logV[m] - X0, y[m], deg=deg) res.append(y - np.polyval(p, logV - X0)) D = np.fmax.reduce(res) if len(res) > 1 else res[0] D = D[np.isfinite(D)] row = dict(variant=name, N=len(D), robust_sigma=mad_sigma(D), n_lt_m03=int((D < -0.301).sum()), n_lt_m06=int((D < -0.602).sum()), n_lt_m10=int((D < -1.0).sum()), frac_lt_m06=float((D < -0.602).mean()), n_gt_p06=int((D > 0.602).sum())) if extra: row.update(extra) return row def main(): d = pd.read_csv(os.path.join(DATA, "alfalfa_all.csv"), low_memory=False) base = (d.HI == 1) & (d.Flag == 1) & d.logMsM.notna() & d.logLi.notna() & d.ba.notna() rows = [] def alf(name, snr=10, dmin=40, bamax=0.70, vmin=80, q0=0.20, bands=("logLi", "logLW1"), deg=2, ew=None, clean=False): cosi2 = np.clip((d.ba ** 2 - q0 ** 2) / (1 - q0 ** 2), 0, 1) V = d.Wc / (2 * np.sqrt(1 - cosi2)) m = base & (d.SNR >= snr) & (d.Dist >= dmin) & (d.ba <= bamax) & (V >= vmin) if ew: m &= d.e_W50 <= ew if clean: m &= d.auto_clean s = d[m] rows.append(stats("ALFALFA-SDSS: " + name, np.log10(V[m]).values, [s[b].values for b in bands], deg=deg)) alf("baseline gold (SNR>=10, D>=40, b/a<=0.7, V>=80, q0=0.2)") alf("catalogue-clean and e_W50<=20 only", ew=20, clean=True) alf("SNR >= 6.5 (ALFALFA code-1 limit)", snr=6.5) alf("SNR >= 15", snr=15) alf("distance >= 60 Mpc", dmin=60) alf("b/a <= 0.5 (inclination > 62 deg)", bamax=0.5) alf("V >= 150 km/s", vmin=150) alf("intrinsic thickness q0 = 0.13", q0=0.13) alf("i band only", bands=("logLi",)) alf("W1 only", bands=("logLW1",)) alf("no internal-extinction correction (i) + W1", bands=("logLi_noint", "logLW1")) alf("linear instead of quadratic relation", deg=1) alf("15-bin running median instead of quadratic relation", deg="rm") alf("15-bin running median, catalogue-clean and e_W50<=20 only", deg="rm", ew=20, clean=True) c = pd.read_csv(os.path.join(DATA, "cf4_resid.csv")) def cf(name, sel, bands=("logLi", "logLW1"), deg=2): s = c[sel] rows.append(stats("CF4: " + name, s.logV.values, [s[b].values for b in bands], deg=deg)) cf("baseline (Vcmb>=3000, V>=80, photometry Q>=3)", np.ones(len(c), bool)) cf("Vcmb >= 5000 km/s", c.Vcmb >= 5000) cf("photometry quality 5 only", (c.QSflag == 5) | (c.QWflag == 5)) cf("inclination error <= 4 deg", c.e_Inc <= 4) cf("V >= 150 km/s", c.V >= 150) cf("i band only", c.logLi.notna(), bands=("logLi",)) cf("W1 only", c.logLW1.notna(), bands=("logLW1",)) cf("only galaxies kept in the CF4 distance table", c.in_table4) cf("linear instead of quadratic relation", np.ones(len(c), bool), deg=1) cf("15-bin running median instead of quadratic relation", np.ones(len(c), bool), deg="rm") u = pd.read_csv(os.path.join(DATA, "union_sample.csv")) for col, lab in (("D", "clipped quadratic (baseline)"), ("D_rm", "15-bin running median")): cl = u[u.auto_clean] rows.append(dict(variant=f"UNION, catalogue-clean: {lab}", N=len(cl), robust_sigma=mad_sigma(cl[col]), n_lt_m03=int((cl[col] < -0.301).sum()), n_lt_m06=int((cl[col] < -0.602).sum()), n_lt_m10=int((cl[col] < -1.0).sum()), frac_lt_m06=float((cl[col] < -0.602).mean()), n_gt_p06=int((cl[col] > 0.602).sum()))) r = pd.DataFrame(rows) r.to_csv(os.path.join(DATA, "robustness.csv"), index=False) pd.set_option("display.width", 250) print(r.round(4).to_string(index=False)) # ---- CF4 membership of ordinary ALFALFA gold galaxies vs apparent magnitude and surface brightness al = pd.read_csv(os.path.join(DATA, "alfalfa_resid.csv"), low_memory=False) g = al[al.gold & al.auto_clean & (al.D.abs() < 0.3) & (al.Vhel < 14000)].copy() # CF4 does not go beyond ~14,000 km/s g["in_cf4"] = g.AGC.isin(c.AGC[c.AGC > 0]) mrows = [] for col, bins, lab in (("imag", [10, 14, 15, 15.5, 16, 16.5, 17, 17.5, 20], "SDSS cmodel i magnitude"), ("mu_i", [17, 20.5, 21, 21.5, 22, 22.5, 23, 23.5, 27], "mean i surface brightness in Petrosian ellipse")): for lo, hi in zip(bins[:-1], bins[1:]): b = g[(g[col] > lo) & (g[col] <= hi)] if len(b) == 0: continue k = int(b.in_cf4.sum()) l68, h68 = clopper_pearson(k, len(b)) mrows.append(dict(quantity=lab, lo=lo, hi=hi, N_gold=len(b), N_in_cf4=k, frac_in_cf4=k / len(b), frac_lo68=l68, frac_hi68=h68)) m = pd.DataFrame(mrows) m.round(3).to_csv(os.path.join(DATA, "cf4_membership.csv"), index=False) print(f" CF4 membership of {len(g)} ordinary ALFALFA gold galaxies (overall {g.in_cf4.mean():.3f}):") print(m.round(3).to_string(index=False)) if __name__ == "__main__": main()