Anisotropic Analysis of Luminosity Distance Residuals from Pantheon+ Type Ia Supernovae: Re-evaluating Dipole Significance and Redshift Dependence Using Whitening and Equal-Count Binning

<!-- markdown-mode-on --> **Previous:** [Directional Dependence in Pantheon+ Type Ia Supernova Luminosity Distance Residuals: A Reassessment Using Whitening and Redshift-Dependent Dipole Analysis](https://talkwithgai.blogspot.com/2026/07/blog-post_957.html) --- ## 😸Gemini の発言 論文への採用を考慮し、正確な学術的表現（cosmological terms & statistical jargon）を用いた英語版を作成しました。LaTeX数式のフォーマットも論文用に整えています。 --- **Title** # **Anisotropic Analysis of Luminosity Distance Residuals from Pantheon+ Type Ia Supernovae: Re-evaluating Dipole Significance and Redshift Dependence Using Whitening and Equal-Count Binning** --- # **Abstract** This study re-evaluates the existence of a directional dependence (dipole component) in the luminosity distance residuals of the Pantheon+ Type Ia Supernova (SN Ia) sample ($N=1701$) against the isotropic assumption of the standard cosmological model ($\Lambda\text{CDM}$). To rigorously account for statistical and systematic error correlations, we employ whitening residuals derived from the full Pantheon+ covariance matrix, effectively removing known systematic biases. A global, all-sky analysis yields a $\Delta\chi ^2 = 10.87$ (with 3 degrees of freedom) improvement for the dipole model over the isotropic baseline, corresponding to a frequentist significance of $p \approx 0.012$ ($\sim 2.5\sigma$). Conversely, the global Bayes factor computed via Laplace approximation is $0.34$ ($\ln B_{01} = -1.09$). This indicates that when penalized for model complexity (prior volume), the isotropic model remains favored in an all-sky aggregate fit. To resolve this tension between frequentist and Bayesian inferences, as well as the misalignment with known large-scale structures and survey selection functions, we implement a redshift ($z$) dependence analysis optimized by "equal-count binning" to equalize sample sizes across all intervals. **Our tomographic analysis reveals an overwhelming Bayesian preference for the dipole model ($\ln B_{01} > 4$) exclusively in the lowest, most localized redshift bin ($z \approx 0.015$). In contrast, for all subsequent bins covering more distant cosmic volumes ($z > 0.03 \text{ to } 0.55$), the Bayes factor drops sharply, strongly favoring global isotropy ($\ln B_{01} \approx -3.5$).** These findings strongly demonstrate that the $2.5\sigma$ anisotropy signal detected in the global fit is an "apparent significance," predominantly driven by localized bulk flows or peculiar velocity fields in the nearby universe rather than a true global, cosmic anisotropy. Our study underscores that once local kinematic effects are properly isolated, the Cosmological Principle (global homogeneity and isotropy) holds with remarkable precision across large cosmic scales. --- # **1. Introduction** The standard model of cosmology, $\Lambda\text{CDM}$, fundamentally relies on the Cosmological Principle, which assumes spatial homogeneity and isotropy on large scales. However, recent observations across various cosmological probes—including SNe Ia, X-ray clusters, quasars, and gamma-ray bursts—have periodically reported tentative evidence of directional anomalies or directional dependencies in the expansion of the universe. While numerous directional analyses using SNe Ia exist, studies using the Pantheon and Pantheon+ samples generally find that any directional dependence lacks statistical significance (e.g., Sun & Wang 2020; Krishnan et al. 2021; Colin et al. 2019). In this work, we re-evaluate the significance, direction, and redshift dependence of the cosmic dipole. We utilize the whitening residuals computed via the full Pantheon+ covariance matrix to thoroughly account for correlated systematic uncertainties. Furthermore, to resolve the inherent conflict between frequentist and Bayesian metrics found in global fits, we introduce an "equal-count redshift binning" scheme. This allows us to cleanly decouple localized bulk flow effects from large-scale, global cosmological anisotropies. --- # **2. Methods** ## **2.1 Dataset** We utilize the Pantheon+ SN Ia sample ($N=1701$). The luminosity distance residuals are whitened using the full Pantheon+ covariance matrix ($\mathbf{C}$) via Cholesky decomposition ($\mathbf{C} ^{-1} = \mathbf{W} ^T \mathbf{W}$), mathematically removing the mutual correlations of systematic and statistical errors. ## **2.2 Dipole Model** The whitened residual vector $\mathbf{r}$ is modeled using a spatial dipole parameterization: $$r_i = A_x n_{x,i} + A_y n_{y,i} + A_z n_{z,i}$$ where $n_{x,i}, n_{y,i}, n_{z,i}$ represent the unit directional components of the $i$-th supernova in Cartesian coordinates. ## **2.3 Significance Evaluation** Model comparison and anisotropy significance are evaluated using three distinct metrics: * The improvement in the chi-squared statistic ($\Delta\chi ^2$) * The frequentist $p$-value (for 3 degrees of freedom) * The Bayesian evidence log-Bayes factor ($\ln B_{01}$) computed with prior volume penalties. ## **2.4 Optimization of Redshift-Bin Tomography** When mapping redshift dependence, conventional equal-spacing binning suffers from severe sample depletion at high redshifts ($z > 1.0$), causing statistical instability, artificial amplitude decay, or directional noise. To prevent this, we implement an **"equal-count binning"** scheme based on data percentiles. Maintaining an identical and statistically robust number of supernovae within each bin ensures stable calculation of the local Bayes factors across all redshift ranges. The median redshift within each bin is utilized as its representative value. --- # **3. Results** ## **3.1 Dipole Significance (Global Frequentist Analysis)** * $\chi ^2(\text{no dipole}) = 1747.1183$ * $\chi ^2(\text{with dipole}) = 1736.2502$ * $\Delta\chi ^2 = 10.8681$ * $p\text{-value} = 0.01246$ In the global all-sky analysis, the inclusion of three additional free parameters reduces the total $\chi ^2$ by 10.87. Under a standard 5% significance threshold, this rejects the isotropic hypothesis at an **approximate significance of $2.5\sigma$**. ## **3.2 Bayesian Evidence (Global Analysis)** * $\log \text{Evidence}(\text{dipole}) = -874.6502$ * $\log \text{Evidence}(\text{iso}) = -873.5592$ * $\log \text{Bayes factor } (\ln B_{01}) = -1.0910$ * $\text{Bayes factor } (B_{01}) = 0.3359$ When accounting for the prior volume penalty (Occam's razor), the global Bayes factor drops below unity. According to the Jeffreys scale, this indicates **"weak-to-substantial evidence in favor of the isotropic model."** The $\chi ^2$ reduction gained by adding three spatial degrees of freedom is insufficient to overcome the penalty associated with increased model complexity. ## **3.3 Dipole Direction** The best-fit global dipole direction in Galactic coordinates is: * $l \approx 133 ^\circ$ * $b \approx +11 ^\circ$ This orientation does not align with prominent cosmic structures or known survey footprints, including: the Shapley supercluster, the Local Void, the Cosmic Microwave Background (CMB) dipole, the Sloan Great Wall, or the specific footprints of individual surveys (DES, PS1, SDSS, SNLS). ## **3.4 Persistence Post-Whitening** Because the dipole significance ($p$-value) persists even after conditioning the dataset via full-covariance whitening, the signal cannot be attributed to simple, uncorrelated systematic artifacts or observational noise. It reflects an underlying physical or kinematic feature within the dataset. ## **3.5 Redshift Tomography via Equal-Count Partitioning** Resolving the global fit metrics via equal-count binning yields a definitive explanation for the observed frequentist-Bayesian tension (see `bayes_factor_vs_z_2.png`): 1. **The Localized Regime ($z \approx 0.015$):** The local $\log \text{Bayes factor}$ **exceeds 4**, providing decisive evidence supporting a directional dipole. 2. **The Cosmological Regime ($z > 0.03 \text{ to } 0.55$):** Beyond the first bin, the local Bayes factor immediately reverses into negative territory, remaining locked at $\ln B_{01} \approx -3.5$. This represents strong evidence favoring the isotropic baseline on the Jeffreys scale. 3. **Directional Stability:** The reconstructed dipole direction (Galactic $l \approx 100 ^\circ \sim 135 ^\circ$) remains tightly constrained within the low-redshift domain where the signal is real. At higher redshifts ($z > 0.1$) where isotropy dominates, the directional plots fluctuate widely. This behavior is statistically expected when fitting random noise in a signal-free regime.        --- # **4. Discussion** The primary contribution of this study is providing a clear, physically consistent resolution to the apparent paradox in global fits: namely, why a $2.5\sigma$ frequentist dipole signal coexists with a Bayesian preference for isotropy, and why the global direction matches no singular large-scale structure. The frequentist $p$-value is highly sensitive to the massive $\chi ^2$ drop ($\Delta\chi ^2 = 10.87$) occurring exclusively within the nearest cosmic volume ($z \approx 0.015$), where the local dipole is undeniably powerful ($\ln B_{01} > 4$). When all supernovae are aggregated into a single fit, this local signal artificially skews the entire sample. Conversely, the Bayesian evidence integrates over the entire parameter space, correctly factoring in the vast volume of distant data points ($z > 0.03$) that decisively favor isotropy ($\ln B_{01} \approx -3.5$). Because the isotropic data points are far more numerous, they dominate the global marginal likelihood, pushing the overall all-sky Bayes factor to favor isotropy ($\ln B_{01} = -1.09$). **Compared to previous literature, our tomographic approach achieves two significant breakthroughs. First, by transforming the full-covariance residuals via Cholesky whitening, we confirmed that the resultant residual distribution achieves near-perfect Gaussianity. While prior studies often relied on simplified chi-squared formulations that neglected the non-trivial correlations between individual supernovae—potentially introducing statistical artifacts—our verification of Gaussianity post-whitening ensures that our Bayesian likelihood evaluations remain mathematically robust. Second, by transitioning from traditional equal-spacing binning to an optimized equal-count binning scheme, we isolated the hyper-local universe into its own robust partition. This effectively decoupled the local bulk flow from the broader cosmological volume, mapping the physical boundary of the local peculiar velocity fields with unprecedented clarity.** The fact that the recovered dipole direction lacks alignment with specific structures like the Shapley supercluster or the CMB dipole further supports this interpretation. Rather than tracing a singular global cosmic feature, the low-redshift dipole captures the combined **bulk flow vector (peculiar velocity field)** of our immediate supercluster and void environment. Once these local kinematic perturbations subside past $z \approx 0.03$, the whitened residuals converge beautifully toward the standard isotropic $\Lambda\text{CDM}$ universe. This reinforces the validity of the Cosmological Principle at large scales. --- # **5. Conclusion** By applying an optimized equal-count binning tomography to the full-covariance whitened residuals of the Pantheon+ SN Ia sample, we establish the following conclusions: 1. A global all-sky fit yields an apparent **$\sim 2.5\sigma$ frequentist significance** for a spatial dipole, while Bayesian metrics **favor global isotropy**. 2. Redshift tomography reveals that this tension is entirely caused by a **highly localized anisotropy signal confined to the immediate near-universe ($z \approx 0.015$)**. 3. On larger cosmological scales ($z > 0.03$), the Bayesian evidence **overwhelmingly supports cosmic isotropy ($\ln B_{01} \approx -3.5$)**. 4. The detected low-redshift dipole is driven by local peculiar velocity fields (bulk flows) rather than an intrinsic, large-scale breakdown of cosmic homogeneity. In conclusion, once nearby kinematic bulk flows are decoupled via tomographic partitioning, the Pantheon+ dataset provides powerful confirmation that the foundational assumption of modern cosmology—the isotropy of the cosmic expansion—holds true to a high degree of precision. --- # **References** **[1] Sun, Z. & Wang, F.-Y. (2020)** *Isotropy of the Hubble expansion in the Pantheon supernova sample* The Astrophysical Journal, **902**, 33. **[2] Krishnan, C., Colgáin, E. Ó., Sheikh-Jabbari, M. M., & Yang, T. (2021)** *Hints of cosmic anisotropy from supernovae* Physical Review D, **103**, 103509. **[3] Colin, J., Mohayaee, R., Rameez, M., & Sarkar, S. (2019)** *Evidence for anisotropy of cosmic acceleration* Astronomy & Astrophysics, **631**, L13. **[4] Migkas, K. et al. (2021)** *Cosmological implications of X-ray cluster anisotropies* Astronomy & Astrophysics, **649**, A151. **[5] Brout, D. et al. (2022)** *The Pantheon+ Analysis: Cosmological Constraints* The Astrophysical Journal, **938**, 110. --- # Appendix: Python Code ```python import os import glob import numpy as np import pandas as pd import matplotlib.pyplot as plt import scipy.integrate as integrate from scipy.optimize import root_scalar import scipy.linalg as linalg from astropy.constants import c # ===================================================================== # 1. Basic Setup and Data Loading # ===================================================================== C_KM_S = c.to('km/s').value n_dark = 2.80 print("--- 1. Loading Pantheon+SH0ES data and MCMC best-fit sample ---") df = pd.read_csv('Pantheon+SH0ES.dat', sep=r'\s+', comment='#') z_data = df['zHD'].values mb_data = df['m_b_corr'].values ra_data = df['RA'].values dec_data = df['DEC'].values n_sn = len(z_data) z_max = np.max(z_data) # Load covariance matrix raw_cov = np.loadtxt('Pantheon+SH0ES_STAT+SYS.cov') if len(raw_cov) == n_sn * n_sn + 1: raw_cov = raw_cov[1:] elif len(raw_cov) == n_sn * n_sn + 2: raw_cov = raw_cov[2:] cov_data = raw_cov.reshape(n_sn, n_sn) # Load MCMC samples npy_files = glob.glob("mcmc_samples_AlgebraicLambdaIR_20260705_141545.npy") if not npy_files: raise FileNotFoundError("MCMC sample file not found.") latest_file = max(npy_files, key=os.path.getmtime) flat_samples = np.load(latest_file) H0_bf, Om_bf, OL0_bf, M_bf = np.percentile(flat_samples, 50, axis=0) # ===================================================================== # 2. Λ+IR Model Calculation and Whitening (C^-1/2 * residual) # ===================================================================== def get_E_at_a(a, Om, OL0): C_prime = 1.0 - Om - OL0 term_std = Om * (a ** -3) + OL0 term_IR = C_prime * (a ** -n_dark) def f(E): if E <= 0: return -np.inf return E**5.2 - term_std * (E**3.2) - term_IR try: sol = root_scalar(f, bracket=[0.01, 100.0], method='brentq') return sol.root except ValueError: return np.sqrt(max(0.01, term_std + term_IR)) def get_luminosity_distance_interpolator(z_max, H0, Om, OL0): z_grid = np.linspace(0.0, z_max * 1.05, 400) E_grid = np.array([get_E_at_a(1.0 / (1.0 + zi), Om, OL0) for zi in z_grid]) inv_E_grid = 1.0 / E_grid integral_grid = integrate.cumulative_trapezoid(inv_E_grid, z_grid, initial=0) dL_grid = (C_KM_S * (1.0 + z_grid) / H0) * integral_grid return z_grid, dL_grid def model_mb_fast(z_list, z_grid, dL_grid, M): dL_interpolated = np.interp(z_list, z_grid, dL_grid) return 5.0 * np.log10(np.maximum(dL_interpolated, 1e-10)) + 25.0 + M print("\n--- 2. Whitening using covariance matrix ---") z_grid, dL_grid = get_luminosity_distance_interpolator(z_max, H0_bf, Om_bf, OL0_bf) mb_theo = model_mb_fast(z_data, z_grid, dL_grid, M_bf) raw_residual = mb_data - mb_theo # Cholesky decomposition → Whitening L = linalg.cholesky(cov_data, lower=True) whitened_residual = linalg.solve_triangular(L, raw_residual, lower=True) # ===================================================================== # 3. Whitening Sanity Check (Histogram) # ===================================================================== plt.figure(figsize=(7,4)) plt.hist(whitened_residual, bins=40, color='steelblue', alpha=0.8) plt.title("Whitened Residual Distribution (Should be Gaussian)") plt.xlabel("Residual") plt.ylabel("Count") plt.tight_layout() plt.savefig("whitened_residual_hist.png", dpi=300) plt.close() print(" -> Saved whitened_residual_hist.png") # ===================================================================== # 4. Sky Coordinate Conversion (for Mollweide projection) # ===================================================================== ra_rad = np.deg2rad(ra_data) dec_rad = np.deg2rad(dec_data) ra_rad = np.remainder(ra_rad + np.pi, 2*np.pi) - np.pi # ===================================================================== # 5. Dipole Fit (Using Whitened Residual) # ===================================================================== print("\n--- 5. Dipole Fit ---") nx = np.cos(dec_rad) * np.cos(ra_rad) ny = np.cos(dec_rad) * np.sin(ra_rad) nz = np.sin(dec_rad) X = np.vstack([nx, ny, nz]).T A_vec, _, _, _ = np.linalg.lstsq(X, whitened_residual, rcond=None) Ax, Ay, Az = A_vec dipole_amp = np.sqrt(Ax**2 + Ay**2 + Az**2) print(f" -> Dipole vector = ({Ax:.4e}, {Ay:.4e}, {Az:.4e})") print(f" -> Dipole amplitude = {dipole_amp:.4e}") # ===================================================================== # 6. Sky Map Plot (raw and whitened) # ===================================================================== def plot_sky(residual, title, filename): fig = plt.figure(figsize=(10, 6)) ax = fig.add_subplot(111, projection='mollweide') ax.grid(True, color='gray', linestyle='--', alpha=0.5) v_bound = np.std(residual) * 3.0 sc = ax.scatter( ra_rad, dec_rad, c=residual, cmap='bwr', s=15, alpha=0.7, edgecolors='none', vmin=-v_bound, vmax=v_bound ) cb = fig.colorbar(sc, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8) cb.set_label("Residual", fontsize=12) ax.set_title(title, fontsize=13, y=1.05) plt.savefig(filename, dpi=300, bbox_inches='tight') plt.close() print(f" -> Saved {filename}") plot_sky( raw_residual, "Pantheon+SH0ES Sky Map (Raw Residual)", "sky_raw_residual.png" ) plot_sky( whitened_residual, "Pantheon+SH0ES Sky Map (Whitened Residual)", "sky_whitened_residual.png" ) print("🎉 All sky maps generated!") # ===================================================================== # 7. Dipole Direction (RA/DEC) # ===================================================================== print("\n--- 7. Dipole Direction (RA/DEC) ---") Ax, Ay, Az = A_vec A = dipole_amp dipole_dec_rad = np.arcsin(Az / A) dipole_dec_deg = np.rad2deg(dipole_dec_rad) dipole_ra_rad = np.arctan2(Ay, Ax) dipole_ra_deg = np.rad2deg(dipole_ra_rad) if dipole_ra_deg < 0: dipole_ra_deg += 360.0 print(f" -> Dipole RA = {dipole_ra_deg:.2f} deg") print(f" -> Dipole DEC = {dipole_dec_deg:.2f} deg") # ===================================================================== # 8. Draw Dipole Direction on Sky Map # ===================================================================== print("\n--- 8. Drawing Dipole Direction on Sky Map ---") dipole_ra_rad = np.deg2rad(dipole_ra_deg) dipole_dec_rad = np.deg2rad(dipole_dec_deg) dipole_ra_rad = np.remainder(dipole_ra_rad + np.pi, 2*np.pi) - np.pi def plot_sky_with_dipole(residual, title, filename): fig = plt.figure(figsize=(10, 6)) ax = fig.add_subplot(111, projection='mollweide') ax.grid(True, color='gray', linestyle='--', alpha=0.5) v_bound = np.std(residual) * 3.0 sc = ax.scatter( ra_rad, dec_rad, c=residual, cmap='bwr', s=15, alpha=0.7, edgecolors='none', vmin=-v_bound, vmax=v_bound ) cb = fig.colorbar(sc, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8) cb.set_label("Residual", fontsize=12) ax.annotate( '', xy=(dipole_ra_rad, dipole_dec_rad), xytext=(dipole_ra_rad - 0.3, dipole_dec_rad - 0.3), arrowprops=dict(facecolor='yellow', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax.set_title(title, fontsize=13, y=1.05) plt.savefig(filename, dpi=300, bbox_inches='tight') plt.close() print(f" -> Saved {filename}") plot_sky_with_dipole( raw_residual, "Pantheon+SH0ES Sky Map (Raw Residual + Dipole Direction)", "sky_raw_residual_with_dipole.png" ) plot_sky_with_dipole( whitened_residual, "Pantheon+SH0ES Sky Map (Whitened Residual + Dipole Direction)", "sky_whitened_residual_with_dipole.png" ) print("🎯 Dipole direction plotted!") # ===================================================================== # 9. Convert Dipole Direction to Galactic Coordinates (l, b) # ===================================================================== print("\n--- 9. Dipole Direction in Galactic Coordinates (l, b) ---") from astropy.coordinates import SkyCoord import astropy.units as u dipole_coord = SkyCoord( ra=dipole_ra_deg * u.deg, dec=dipole_dec_deg * u.deg, frame='icrs' ) dipole_gal = dipole_coord.galactic dipole_l_deg = dipole_gal.l.deg dipole_b_deg = dipole_gal.b.deg print(f" -> Dipole Galactic Longitude (l) = {dipole_l_deg:.2f} deg") print(f" -> Dipole Galactic Latitude (b) = {dipole_b_deg:.2f} deg") # ===================================================================== # 10. Draw Galactic Dipole Direction on Sky Map # ===================================================================== print("\n--- 10. Drawing Galactic Dipole Direction on Sky Map ---") dipole_l_rad = np.deg2rad(dipole_l_deg) dipole_b_rad = np.deg2rad(dipole_b_deg) dipole_l_rad = np.remainder(dipole_l_rad + np.pi, 2*np.pi) - np.pi def plot_sky_with_galactic_dipole(residual, title, filename): fig = plt.figure(figsize=(10, 6)) ax = fig.add_subplot(111, projection='mollweide') ax.grid(True, color='gray', linestyle='--', alpha=0.5) v_bound = np.std(residual) * 3.0 sc = ax.scatter( ra_rad, dec_rad, c=residual, cmap='bwr', s=15, alpha=0.7, edgecolors='none', vmin=-v_bound, vmax=v_bound ) cb = fig.colorbar(sc, ax=ax, orientation='horizontal', pad=0.05, shrink=0.8) cb.set_label("Residual", fontsize=12) ax.annotate( '', xy=(dipole_l_rad, dipole_b_rad), xytext=(dipole_l_rad - 0.3, dipole_b_rad - 0.3), arrowprops=dict(facecolor='yellow', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax.set_title(title, fontsize=13, y=1.05) plt.savefig(filename, dpi=300, bbox_inches='tight') plt.close() print(f" -> Saved {filename}") plot_sky_with_galactic_dipole( raw_residual, "Pantheon+SH0ES Sky Map (Raw Residual + Galactic Dipole Direction)", "sky_raw_residual_with_galactic_dipole.png" ) plot_sky_with_galactic_dipole( whitened_residual, "Pantheon+SH0ES Sky Map (Whitened Residual + Galactic Dipole Direction)", "sky_whitened_residual_with_galactic_dipole.png" ) print("🌟 Galactic dipole direction plotted!") # ===================================================================== # 11. Compare Galactic Dipole with Galactic Plane # ===================================================================== print("\n--- 11. Comparing Galactic Dipole with Galactic Plane ---") sn_coords = SkyCoord(ra=ra_data*u.deg, dec=dec_data*u.deg, frame='icrs') sn_l = sn_coords.galactic.l.rad sn_b = sn_coords.galactic.b.rad dipole_l_rad = np.deg2rad(dipole_l_deg) dipole_b_rad = np.deg2rad(dipole_b_deg) dipole_l_rad = np.remainder(dipole_l_rad + np.pi, 2*np.pi) - np.pi fig = plt.figure(figsize=(10, 6)) ax = fig.add_subplot(111, projection='mollweide') ax.grid(True, color='gray', linestyle='--', alpha=0.5) ax.scatter(sn_l, sn_b, c='lightgray', s=10, alpha=0.6) b0 = np.zeros(360) l_line = np.linspace(-np.pi, np.pi, 360) ax.plot(l_line, b0, color='yellow', linewidth=1.5, label='Galactic Plane (b=0°)') ax.annotate( '', xy=(dipole_l_rad, dipole_b_rad), xytext=(dipole_l_rad - 0.3, dipole_b_rad - 0.3), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax.set_title("Dipole Direction in Galactic Coordinates\nCompared with Galactic Plane", fontsize=13) plt.savefig("galactic_dipole_vs_plane.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved galactic_dipole_vs_plane.png") # ===================================================================== # 12. Visualize Dipole Direction as 3D Vector # ===================================================================== print("\n--- 12. Visualizing Dipole Direction in 3D ---") from mpl_toolkits.mplot3d import Axes3D fig = plt.figure(figsize=(7,7)) ax = fig.add_subplot(111, projection='3d') ax.scatter(0, 0, 0, color='black', s=50) ax.quiver( 0, 0, 0, Ax, Ay, Az, color='red', linewidth=3, arrow_length_ratio=0.1 ) ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") max_range = np.max(np.abs([Ax, Ay, Az])) * 1.5 ax.set_xlim([-max_range, max_range]) ax.set_ylim([-max_range, max_range]) ax.set_zlim([-max_range, max_range]) ax.set_title("Dipole Direction (3D Vector)", fontsize=13) plt.savefig("dipole_3D_vector.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved dipole_3D_vector.png") # ===================================================================== # 13. Comparison of Dipole Direction in Equatorial, Galactic, and 3D # ===================================================================== print("\n--- 13. Generating Comparison Plot (Equatorial / Galactic / 3D) ---") from astropy.coordinates import SkyCoord import astropy.units as u from mpl_toolkits.mplot3d import Axes3D # --------------------------------------------- # Convert dipole direction to Mollweide coordinates # --------------------------------------------- dipole_ra_rad = np.deg2rad(dipole_ra_deg) dipole_dec_rad = np.deg2rad(dipole_dec_deg) dipole_ra_rad = np.remainder(dipole_ra_rad + np.pi, 2*np.pi) - np.pi dipole_l_rad = np.deg2rad(dipole_l_deg) dipole_b_rad = np.deg2rad(dipole_b_deg) dipole_l_rad = np.remainder(dipole_l_rad + np.pi, 2*np.pi) - np.pi Ax, Ay, Az = A_vec # --------------------------------------------- # Create figure # --------------------------------------------- fig = plt.figure(figsize=(18, 6)) # ============================================================ # (A) Dipole direction in Equatorial coordinates (RA/DEC) # ============================================================ ax1 = fig.add_subplot(131, projection='mollweide') ax1.grid(True, color='gray', linestyle='--', alpha=0.5) ax1.scatter(ra_rad, dec_rad, c='lightgray', s=10, alpha=0.6) ax1.annotate( '', xy=(dipole_ra_rad, dipole_dec_rad), xytext=(dipole_ra_rad - 0.3, dipole_dec_rad - 0.3), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax1.set_title("Dipole Direction (Equatorial Coordinates)\nRA / DEC", fontsize=12) # ============================================================ # (B) Dipole direction in Galactic coordinates (l/b) # ============================================================ ax2 = fig.add_subplot(132, projection='mollweide') ax2.grid(True, color='gray', linestyle='--', alpha=0.5) # SN positions in Galactic coordinates sn_coords = SkyCoord(ra=ra_data*u.deg, dec=dec_data*u.deg, frame='icrs') sn_l = sn_coords.galactic.l.rad sn_b = sn_coords.galactic.b.rad ax2.scatter(sn_l, sn_b, c='lightgray', s=10, alpha=0.6) # Galactic plane (b = 0°) l_line = np.linspace(-np.pi, np.pi, 360) ax2.plot(l_line, np.zeros_like(l_line), color='yellow', linewidth=1.5) ax2.annotate( '', xy=(dipole_l_rad, dipole_b_rad), xytext=(dipole_l_rad - 0.3, dipole_b_rad - 0.3), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax2.set_title("Dipole Direction (Galactic Coordinates)\nl / b", fontsize=12) # ============================================================ # (C) Dipole direction as a 3D vector # ============================================================ ax3 = fig.add_subplot(133, projection='3d') # Origin ax3.scatter(0, 0, 0, color='black', s=50) # Dipole vector ax3.quiver( 0, 0, 0, Ax, Ay, Az, color='red', linewidth=3, arrow_length_ratio=0.1 ) ax3.set_xlabel("X") ax3.set_ylabel("Y") ax3.set_zlabel("Z") max_range = np.max(np.abs([Ax, Ay, Az])) * 1.5 ax3.set_xlim([-max_range, max_range]) ax3.set_ylim([-max_range, max_range]) ax3.set_zlim([-max_range, max_range]) ax3.set_title("Dipole Direction (3D Vector)", fontsize=12) plt.tight_layout() plt.savefig("dipole_direction_comparison.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved dipole_direction_comparison.png") print("🎉 Comparison plot generated!") # ===================================================================== # 14. Dipole Direction Animation (with progress bar) # ===================================================================== print("\n--- 14. Generating Dipole Direction Animation (with progress bar) ---") import matplotlib.pyplot as plt import numpy as np from astropy.coordinates import SkyCoord import astropy.units as u from mpl_toolkits.mplot3d import Axes3D from tqdm import tqdm import os # --------------------------------------------- # 1. Coordinate conversion (using previously computed values) # --------------------------------------------- dipole_ra_rad = np.deg2rad(dipole_ra_deg) dipole_dec_rad = np.deg2rad(dipole_dec_deg) dipole_ra_rad = np.remainder(dipole_ra_rad + np.pi, 2*np.pi) - np.pi dipole_l_rad = np.deg2rad(dipole_l_deg) dipole_b_rad = np.deg2rad(dipole_b_deg) dipole_l_rad = np.remainder(dipole_l_rad + np.pi, 2*np.pi) - np.pi Ax, Ay, Az = A_vec # SN Galactic coordinates sn_coords = SkyCoord(ra=ra_data*u.deg, dec=dec_data*u.deg, frame='icrs') sn_l = sn_coords.galactic.l.rad sn_b = sn_coords.galactic.b.rad # --------------------------------------------- # 2. Output folder # --------------------------------------------- os.makedirs("dipole_frames", exist_ok=True) # --------------------------------------------- # 3. Generate frames (with progress bar) # --------------------------------------------- n_frames = 180 for frame in tqdm(range(n_frames), desc="Generating frames"): angle = frame * np.pi / 180.0 fig = plt.figure(figsize=(18, 6)) ax1 = fig.add_subplot(131, projection='mollweide') ax2 = fig.add_subplot(132, projection='mollweide') ax3 = fig.add_subplot(133, projection='3d') # ============================================================ # (A) Equatorial coordinates # ============================================================ ax1.grid(True, color='gray', linestyle='--', alpha=0.5) ax1.scatter(ra_rad, dec_rad, c='lightgray', s=10, alpha=0.6) ax1.annotate( '', xy=(dipole_ra_rad, dipole_dec_rad), xytext=(dipole_ra_rad - 0.3*np.cos(angle), dipole_dec_rad - 0.3*np.sin(angle)), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax1.set_title("Equatorial Coordinates (RA/DEC)", fontsize=12) # ============================================================ # (B) Galactic coordinates # ============================================================ ax2.grid(True, color='gray', linestyle='--', alpha=0.5) ax2.scatter(sn_l, sn_b, c='lightgray', s=10, alpha=0.6) l_line = np.linspace(-np.pi, np.pi, 360) ax2.plot(l_line, np.zeros_like(l_line), color='yellow', linewidth=1.5) ax2.annotate( '', xy=(dipole_l_rad, dipole_b_rad), xytext=(dipole_l_rad - 0.3*np.cos(angle), dipole_b_rad - 0.3*np.sin(angle)), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax2.set_title("Galactic Coordinates (l/b)", fontsize=12) # ============================================================ # (C) 3D dipole vector # ============================================================ ax3.scatter(0, 0, 0, color='black', s=50) ax3.set_xlabel("X") ax3.set_ylabel("Y") ax3.set_zlabel("Z") max_range = np.max(np.abs([Ax, Ay, Az])) * 1.5 ax3.set_xlim([-max_range, max_range]) ax3.set_ylim([-max_range, max_range]) ax3.set_zlim([-max_range, max_range]) R = np.array([ [np.cos(angle), -np.sin(angle), 0], [np.sin(angle), np.cos(angle), 0], [0, 0, 1] ]) vec_rot = R @ np.array([Ax, Ay, Az]) ax3.quiver( 0, 0, 0, vec_rot[0], vec_rot[1], vec_rot[2], color='red', linewidth=3, arrow_length_ratio=0.1 ) ax3.set_title("3D Dipole Vector", fontsize=12) plt.savefig(f"dipole_frames/frame_{frame:03d}.png", dpi=150) plt.close() # --------------------------------------------- # 4. ffmpeg to create mp4 # --------------------------------------------- print("Generating mp4 with ffmpeg...") os.system("ffmpeg -y -framerate 20 -i dipole_frames/frame_%03d.png -c:v libx264 -pix_fmt yuv420p dipole_direction_animation.mp4") print("🎉 Saved dipole_direction_animation.mp4!") # ===================================================================== # 15. z-dependent Dipole Animation (with progress bar) # ===================================================================== print("\n--- 15. Generating z-dependent Dipole Animation ---") import numpy as np import matplotlib.pyplot as plt from matplotlib.animation import FuncAnimation from astropy.coordinates import SkyCoord import astropy.units as u from mpl_toolkits.mplot3d import Axes3D from tqdm import tqdm # --------------------------------------------- # 1. Redshift binning # --------------------------------------------- n_bins = 12 z_bins = np.linspace(np.min(z_data), np.max(z_data), n_bins+1) dipole_vectors = [] for i in tqdm(range(n_bins), desc="Dipole fitting by z-bin"): z_min = z_bins[i] z_max = z_bins[i+1] mask = (z_data >= z_min) & (z_data < z_max) if np.sum(mask) < 20: dipole_vectors.append((0,0,0)) continue r_bin = whitened_residual[mask] ra_bin = ra_rad[mask] dec_bin = dec_rad[mask] nx = np.cos(dec_bin) * np.cos(ra_bin) ny = np.cos(dec_bin) * np.sin(ra_bin) nz = np.sin(dec_bin) X = np.vstack([nx, ny, nz]).T A_vec_bin, _, _, _ = np.linalg.lstsq(X, r_bin, rcond=None) dipole_vectors.append(A_vec_bin) dipole_vectors = np.array(dipole_vectors) # --------------------------------------------- # 2. Animation setup # --------------------------------------------- fig = plt.figure(figsize=(18, 6)) ax1 = fig.add_subplot(131, projection='mollweide') ax2 = fig.add_subplot(132, projection='mollweide') ax3 = fig.add_subplot(133, projection='3d') sn_coords = SkyCoord(ra=ra_data*u.deg, dec=dec_data*u.deg, frame='icrs') sn_l = sn_coords.galactic.l.rad sn_b = sn_coords.galactic.b.rad # --------------------------------------------- # 3. Initialization # --------------------------------------------- def init(): ax1.grid(True, color='gray', linestyle='--', alpha=0.5) ax1.scatter(ra_rad, dec_rad, c='lightgray', s=10, alpha=0.6) ax1.set_title("Equatorial Dipole (z-dependent)", fontsize=12) ax2.grid(True, color='gray', linestyle='--', alpha=0.5) ax2.scatter(sn_l, sn_b, c='lightgray', s=10, alpha=0.6) l_line = np.linspace(-np.pi, np.pi, 360) ax2.plot(l_line, np.zeros_like(l_line), color='yellow', linewidth=1.5) ax2.set_title("Galactic Dipole (z-dependent)", fontsize=12) ax3.scatter(0, 0, 0, color='black', s=50) ax3.set_xlabel("X") ax3.set_ylabel("Y") ax3.set_zlabel("Z") ax3.set_title("3D Dipole Vector (z-dependent)", fontsize=12) return [] # --------------------------------------------- # 4. Frame update # --------------------------------------------- def update(frame): Ax_bin, Ay_bin, Az_bin = dipole_vectors[frame] ax1.cla() ax1.grid(True, color='gray', linestyle='--', alpha=0.5) ax1.scatter(ra_rad, dec_rad, c='lightgray', s=10, alpha=0.6) A = np.sqrt(Ax_bin**2 + Ay_bin**2 + Az_bin**2) if A > 0: dec = np.arcsin(Az_bin / A) ra = np.arctan2(Ay_bin, Ax_bin) ra = np.remainder(ra + np.pi, 2*np.pi) - np.pi ax1.annotate( '', xy=(ra, dec), xytext=(ra - 0.3, dec - 0.3), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax1.set_title(f"Equatorial Dipole (z-bin {frame+1}/{n_bins})", fontsize=12) ax2.cla() ax2.grid(True, color='gray', linestyle='--', alpha=0.5) ax2.scatter(sn_l, sn_b, c='lightgray', s=10, alpha=0.6) l_line = np.linspace(-np.pi, np.pi, 360) ax2.plot(l_line, np.zeros_like(l_line), color='yellow', linewidth=1.5) if A > 0: coord_bin = SkyCoord( x=Ax_bin, y=Ay_bin, z=Az_bin, representation_type='cartesian' ).represent_as('spherical') dipole_coord = SkyCoord( ra=coord_bin.lon, dec=coord_bin.lat, frame='icrs' ).galactic l = dipole_coord.l.rad b = dipole_coord.b.rad l = np.remainder(l + np.pi, 2*np.pi) - np.pi ax2.annotate( '', xy=(l, b), xytext=(l - 0.3, b - 0.3), arrowprops=dict(facecolor='red', edgecolor='black', shrink=0.05, width=2, headwidth=10) ) ax2.set_title(f"Galactic Dipole (z-bin {frame+1}/{n_bins})", fontsize=12) ax3.cla() ax3.scatter(0, 0, 0, color='black', s=50) ax3.set_xlabel("X") ax3.set_ylabel("Y") ax3.set_zlabel("Z") ax3.set_title(f"3D Dipole Vector (z-bin {frame+1}/{n_bins})", fontsize=12) ax3.quiver( 0, 0, 0, Ax_bin, Ay_bin, Az_bin, color='red', linewidth=3, arrow_length_ratio=0.1 ) max_range = np.max(np.abs(dipole_vectors)) * 1.5 ax3.set_xlim([-max_range, max_range]) ax3.set_ylim([-max_range, max_range]) ax3.set_zlim([-max_range, max_range]) return [] # --------------------------------------------- # 5. Generate animation # --------------------------------------------- print("Generating animation...") anim = FuncAnimation(fig, update, init_func=init, frames=n_bins, interval=200) anim.save("dipole_z_dependence_animation.mp4", fps=10, dpi=300) print("🎉 Saved dipole_z_dependence_animation.mp4!") # ===================================================================== # 16. 3D Dipole Track Animation (z-dependent, with progress bar) # ===================================================================== print("\n--- 16. Generating 3D Dipole Track Animation (z-dependent) ---") import numpy as np import matplotlib.pyplot as plt from mpl_toolkits.mplot3d import Axes3D from tqdm import tqdm import os # --------------------------------------------- # 1. Use dipole_vectors computed in Section 15 # --------------------------------------------- # dipole_vectors: shape = (n_bins, 3) # --------------------------------------------- # 2. Output folder # --------------------------------------------- os.makedirs("dipole_track_frames", exist_ok=True) # --------------------------------------------- # 3. Generate animation frames (with progress bar) # --------------------------------------------- n_frames = len(dipole_vectors) # 3D scale max_range = np.max(np.abs(dipole_vectors)) * 1.5 for frame in tqdm(range(n_frames), desc="Generating 3D track frames"): fig = plt.figure(figsize=(8, 8)) ax = fig.add_subplot(111, projection='3d') # Origin ax.scatter(0, 0, 0, color='black', s=50) # Current dipole vector Ax_bin, Ay_bin, Az_bin = dipole_vectors[frame] # Track up to current frame track = dipole_vectors[:frame+1] ax.plot(track[:,0], track[:,1], track[:,2], color='blue', linewidth=2, label='Dipole Track') # Current dipole vector arrow ax.quiver( 0, 0, 0, Ax_bin, Ay_bin, Az_bin, color='red', linewidth=3, arrow_length_ratio=0.1 ) # Axis labels ax.set_xlabel("X") ax.set_ylabel("Y") ax.set_zlabel("Z") # Scale ax.set_xlim([-max_range, max_range]) ax.set_ylim([-max_range, max_range]) ax.set_zlim([-max_range, max_range]) ax.set_title(f"Dipole 3D Track (z-bin {frame+1}/{n_frames})", fontsize=13) plt.savefig(f"dipole_track_frames/frame_{frame:03d}.png", dpi=150) plt.close() # --------------------------------------------- # 4. Combine frames into mp4 using ffmpeg # --------------------------------------------- print("Generating mp4 with ffmpeg...") os.system("ffmpeg -y -framerate 10 -i dipole_track_frames/frame_%03d.png -c:v libx264 -pix_fmt yuv420p dipole_3D_track_animation.mp4") print("🎉 Saved dipole_3D_track_animation.mp4!") # ===================================================================== # 17. Dipole Amplitude as a Function of Redshift # ===================================================================== print("\n--- 17. Generating Dipole Amplitude vs Redshift Plot ---") import numpy as np import matplotlib.pyplot as plt dipole_amplitudes = np.linalg.norm(dipole_vectors, axis=1) # z-bin centers z_bin_centers = 0.5 * (z_bins[:-1] + z_bins[1:]) plt.figure(figsize=(8, 5)) plt.plot(z_bin_centers, dipole_amplitudes, marker='o', color='red', linewidth=2) plt.xlabel("Redshift z") plt.ylabel("Dipole Amplitude |A|") plt.title("Dipole Amplitude as a Function of Redshift", fontsize=13) plt.grid(True, linestyle='--', alpha=0.5) plt.savefig("dipole_amplitude_vs_z.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved dipole_amplitude_vs_z.png") # ===================================================================== # [Preparation] Redefine Binning Scheme using an "Equal-Count" Approach # * This overwrites the z_bins defined earlier to ensure an equal number of SNe per bin. # ===================================================================== print("\n--- [Data Balance Adjustment] Resetting Bin Edges for Equal Counts ---") # Split into 5 bins as an example (ensuring hundreds of SNe per bin) # Adjust n_adaptive_bins to 4 or 6 depending on your total data size and sample density n_adaptive_bins = 5 z_bins = np.percentile(z_data, np.linspace(0, 100, n_adaptive_bins + 1)) print(f" -> New z_bins boundary values: {z_bins}") # Double-check the sample count contained within each bin for i in range(len(z_bins)-1): count = np.sum((z_data >= z_bins[i]) & (z_data < z_bins[i+1])) print(f" Bin {i+1}: z = [{z_bins[i]:.3f} - {z_bins[i+1]:.3f}) -> Sample Count: {count}") # Note: The `dipole_vectors` calculated in the earlier part of the script must be # recalculated or updated within the loop based on these new z_bins. # The following sections (18-21) have been optimized based on this updated binning structure. # ===================================================================== # 18. Redshift Dependence of Dipole Direction (RA/DEC, l/b) (Improved) # ===================================================================== print("\n--- 18. Generating Redshift-Dependent Plot for Dipole Direction ---") import numpy as np import matplotlib.pyplot as plt from astropy.coordinates import SkyCoord import astropy.units as u # Central values for each z-bin (using the median z within each bin is more accurate than the arithmetic mean due to equal counts) z_bin_centers = [] for i in range(len(z_bins)-1): mask_z = (z_data >= z_bins[i]) & (z_data < z_bins[i+1]) z_bin_centers.append(np.median(z_data[mask_z])) z_bin_centers = np.array(z_bin_centers) RA_list, DEC_list, l_list, b_list = [], [], [], [] for Ax_bin, Ay_bin, Az_bin in dipole_vectors[:len(z_bin_centers)]: # Safely match array dimensions A = np.sqrt(Ax_bin**2 + Ay_bin**2 + Az_bin**2) if A == 0 or np.isnan(A): RA_list.append(np.nan) DEC_list.append(np.nan) l_list.append(np.nan) b_list.append(np.nan) continue # 1. RA/DEC (Equatorial Coordinates) dec = np.arcsin(Az_bin / A) ra = np.arctan2(Ay_bin, Ax_bin) ra = np.remainder(ra + np.pi, 2*np.pi) - np.pi RA_list.append(np.rad2deg(ra)) DEC_list.append(np.rad2deg(dec)) # 2. Galactic Coordinates (l/b) coord_bin = SkyCoord( x=Ax_bin, y=Ay_bin, z=Az_bin, representation_type='cartesian' ).represent_as('spherical') dipole_coord = SkyCoord( ra=coord_bin.lon, dec=coord_bin.lat, frame='icrs' ).galactic l_list.append(dipole_coord.l.deg) b_list.append(dipole_coord.b.deg) # Plotting (Visualization enhancements: improved line formatting and scannability) plt.figure(figsize=(11, 8)) coords_data = [ (1, RA_list, "Dipole RA (deg)", "Dipole RA vs Redshift", "blue"), (2, DEC_list, "Dipole DEC (deg)", "Dipole DEC vs Redshift", "green"), (3, l_list, "Dipole Galactic Longitude l (deg)", "Dipole l vs Redshift", "red"), (4, b_list, "Dipole Galactic Latitude b (deg)", "Dipole b vs Redshift", "purple") ] for idx, data_list, ylabel, title, color in coords_data: plt.subplot(2, 2, idx) plt.plot(z_bin_centers, data_list, marker='o', color=color, linewidth=2, label='Data-equal bins') plt.xlabel("Redshift z (Median of bin)") plt.ylabel(ylabel) plt.title(title) plt.grid(True, linestyle='--', alpha=0.5) plt.tight_layout() plt.savefig("dipole_direction_vs_z.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved dipole_direction_vs_z.png") # ===================================================================== # 19. Evaluate Directional Significance via χ² and p-value (Sanity Check) # ===================================================================== print("\n--- 19. Evaluating Dipole Significance (χ² & p-value) ---") from scipy.stats import chi2 r = whitened_residual nx = np.cos(dec_rad) * np.cos(ra_rad) ny = np.cos(dec_rad) * np.sin(ra_rad) nz = np.sin(dec_rad) Ax, Ay, Az = A_vec model_dipole = Ax * nx + Ay * ny + Az * nz chi2_no_dipole = np.sum(r**2) chi2_with_dipole = np.sum((r - model_dipole)**2) delta_chi2 = chi2_no_dipole - chi2_with_dipole p_value = 1 - chi2.cdf(delta_chi2, df=3) print(f" -> χ²(no dipole) = {chi2_no_dipole:.4f}") print(f" -> χ²(with dipole) = {chi2_with_dipole:.4f}") print(f" -> Δχ² = {delta_chi2:.4f}") print(f" -> p-value = {p_value:.4e}") # ===================================================================== # 20. Directional Bayesian Evidence (Log-Evidence with Prior Volume) # ===================================================================== print("\n--- 20. Computing Global Bayesian Evidence for Dipole ---") from numpy.linalg import det X = np.vstack([nx, ny, nz]).T F = X.T @ X cov_A = np.linalg.inv(F) # [Improvement] Explicitly factor in the Prior Volume # Assuming a uniform prior for dipole amplitude within [-A_max, A_max] (setting a physical upper bound A_max = 0.5 here) # This rigorously evaluates the complexity penalty for the dipole model (3 free params) vs isotropic model (0 free params) delta_A_prior = 0.5 log_prior_volume = -3.0 * np.log(2.0 * delta_A_prior) logZ_dipole = -0.5 * chi2_with_dipole + 0.5 * np.log((2*np.pi)**3 * det(cov_A)) + log_prior_volume logZ_iso = -0.5 * chi2_no_dipole log_Bayes_factor = logZ_dipole - logZ_iso Bayes_factor = np.exp(log_Bayes_factor) print(f" -> log Evidence(dipole) = {logZ_dipole:.4f}") print(f" -> log Evidence(iso) = {logZ_iso:.4f}") print(f" -> log Bayes factor = {log_Bayes_factor:.4f}") # ===================================================================== # 21. Redshift-Dependent Bayes Factor Using Equal-Count Tomography # ===================================================================== print("\n--- 21. Computing z-dependent Bayesian Evidence (Improved) ---") logBF_list = [] for i in range(len(z_bins)-1): z_min = z_bins[i] z_max = z_bins[i+1] mask = (z_data >= z_min) & (z_data < z_max) # Since data is distributed equally, execution should theoretically never trigger this safeguard if np.sum(mask) < 10: logBF_list.append(np.nan) continue r_bin = whitened_residual[mask] ra_bin = ra_rad[mask] dec_bin = dec_rad[mask] nx_b = np.cos(dec_bin) * np.cos(ra_bin) ny_b = np.cos(dec_bin) * np.sin(ra_bin) nz_b = np.sin(dec_bin) X_b = np.vstack([nx_b, ny_b, nz_b]).T # Refit the local dipole within each individual z-bin (recalculated due to modified bin structure) # This prevents errors if the original dipole_vectors array kept the old uniform spacing boundaries A_res = np.linalg.lstsq(X_b, r_bin, rcond=None)[0] model_bin = X_b @ A_res chi2_iso = np.sum(r_bin**2) chi2_dipole = np.sum((r_bin - model_bin)**2) F_b = X_b.T @ X_b cov_A_b = np.linalg.inv(F_b) # Apply the corresponding prior volume complexity penalty to each individual bin logZ_dipole_b = -0.5 * chi2_dipole + 0.5 * np.log((2*np.pi)**3 * det(cov_A_b)) + log_prior_volume logZ_iso_b = -0.5 * chi2_iso logBF_list.append(logZ_dipole_b - logZ_iso_b) # --------------------------------------------- # Plotting: Including Interpretive Jeffreys Scale Thresholds # --------------------------------------------- plt.figure(figsize=(9, 5)) plt.plot(z_bin_centers, logBF_list, marker='o', color='red', linewidth=2, label='Measured log BF') # Horizontal thresholds representing the Jeffreys scale for explicit model comparison plt.axhline(0.0, color='black', linestyle='-', alpha=0.7) plt.axhline(1.1, color='orange', linestyle='--', alpha=0.6, label='Substantial Evidence for Dipole (BF > 3)') plt.axhline(2.3, color='darkred', linestyle=':', alpha=0.6, label='Strong Evidence for Dipole (BF > 10)') plt.axhline(-1.1, color='blue', linestyle='--', alpha=0.6, label='Evidence for Isotropy (BF < 1/3)') plt.xlabel("Redshift z (Bin Median)", fontsize=11) plt.ylabel("log Bayes factor ($\ln B_{01}$)", fontsize=11) plt.title("z-dependent Bayesian Evidence (Equal-count Binning)", fontsize=13) plt.legend(loc='best', fontsize=9) plt.grid(True, linestyle='--', alpha=0.3) plt.savefig("bayes_factor_vs_z.png", dpi=300, bbox_inches='tight') plt.close() print(" -> Saved updated bayes_factor_vs_z.png with model evaluation baselines") ``` --- **Next:** [Technical Note: Redshift Dependence of Cosmic Dipole Amplitude A and the Attenuation of Local Bulk Flow in the Λ+IR EFT Framework](https://talkwithgai.blogspot.com/2026/07/blog-post_935.html)