astronomy-inference

# Astronomy inference Astronomy overloads its conventions: magnitudes are logarithmic and zero-point dependent, times have scales as well as formats, coordinates have frames and epochs, and a periodogram peak is a candidate rather than a detection. Fix the conventions first, then infer, then report the priors and diagnostics alongside the posterior. ## astropy foundations 1. Attach units at input and keep them until the final print: `astropy.units` Quantities with `.to(u.Jy)`, spectral conversions via `equivalencies=u.spectral()` and flux density conversions via `u.spectral_density(wavelength)`; constants from `astropy.constants`. Stripping units early is how a factor of 1e-23 disappears. 2. Coordinates: `SkyCoord(ra, dec, unit="deg", frame="icrs")`; ICRS, FK5 (J2000) and Galactic are different frames. Angular separation is `c1.separation(c2)`, never a Pythagorean difference in RA and Dec (RA shrinks by cos Dec). Cross-match with `match_to_catalog_sky` and report the match radius. High proper-motion targets need `obstime` and `apply_space_motion` to the observation epoch. 3. Time: `Time("2023-01-01T00:00:00", format="isot", scale="utc")`; scales UTC, TAI, TT and TDB differ by tens of seconds, formats `jd`, `mjd`, `isot` are only representations (MJD = JD - 2400000.5). Timing work uses BJD_TDB: `ltt = t.light_travel_time(coord, kind="barycentric", location=site)` then `(t.tdb + ltt).jd`. Heliocentric JD differs from barycentric by several seconds; UTC differs from TDB by more than a minute. State the scale and reference frame of every time column. 4. FITS: `with fits.open(path) as hdul: hdul.info()`, then `hdul[i].header` and `hdul[i].data`; read `BUNIT`, `BSCALE`, `BZERO`, `EXPTIME`, `MJD-OBS`. Tables load with `astropy.table.Table.read(path, hdu=1)`. FITS data are big-endian; convert before handing arrays to libraries that assume native order. WCS: `w = WCS(header)`, `w.pixel_to_world(x, y)` and `w.world_to_pixel(coord)` use zero-based pixels, while FITS headers are one-based (`origin` argument in `wcs_pix2world`). Pixel scale via `astropy.wcs.utils.proj_plane_pixel_scales(w)`. ## Photometry 5. Magnitudes are m = -2.5 log10(F / F0) and the zero point F0 defines the system: AB uses 3631 Jy (`flux.to(u.ABmag)` for a spectral flux density), Vega uses the spectrum of Vega per filter (offsets from AB range from about 0 in V to about 1.85 in Ks), ST uses a flat F_lambda (`u.STmag`). Mixing systems is a 0.1 to 2 magnitude error. Average fluxes, not magnitudes; convert errors with sigma_m = 1.0857 sigma_F / F. Distance modulus mu = m - M = 5 log10(d / 10 pc). Report system, filter, zero point, airmass and extinction corrections, and aperture. 6. Period finding: `ls = astropy.timeseries.LombScargle(t, y, dy)` then `freq, power = ls.autopower(samples_per_peak=10)` with `minimum_frequency` and `maximum_frequency` chosen from the cadence and baseline; the grid must resolve peak widths of order 1 / baseline, so oversample. Significance via `ls.false_alarm_probability(power.max(), method="baluev")` or `method="bootstrap"`; then check aliases (1 cycle per day for ground-based data, harmonics at P/2 and 2P), phase-fold at the candidate period and overlay `ls.model(t_fit, best_frequency)`. Period uncertainty comes from the peak width or a bootstrap, not from the grid spacing. `BoxLeastSquares` is the transit analogue. ## Bayesian inference 7. Write down the likelihood (Gaussian with per-point sigma, plus a jitter term `ln s` when residuals exceed the errors), the prior for every parameter with bounds and shape (uniform versus log-uniform for scale parameters), and the parameterization (sample `sqrt(e) cos w, sqrt(e) sin w` rather than `e, w`). These go in the report verbatim. 8. emcee: `sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=(x, y, yerr))`, `sampler.run_mcmc(p0, nsteps, progress=True)`, then `tau = sampler.get_autocorr_time(tol=0)`. A chain is usable when its length exceeds about 50 tau for every parameter (emcee's default `tol=50` enforces this); discard a burn-in of a few tau and thin by about tau / 2 with `sampler.get_chain(discard=..., thin=..., flat=True)`. Check `sampler.acceptance_fraction` (roughly 0.2 to 0.5). Walkers in one ensemble are not independent chains, so R-hat needs several independent runs; compute it with `arviz.rhat` and effective sample size with `arviz.ess`. Report tau, N / tau and R-hat. 9. dynesty: `sampler = dynesty.NestedSampler(loglike, prior_transform, ndim, nlive=500)`, `sampler.run_nested(dlogz=0.01, checkpoint_file="run.save")`, `res = sampler.results`; evidence `res.logz[-1]` with `res.logzerr[-1]`, equal-weight posterior draws from `res.samples_equal()` (or `dynesty.utils.resample_equal` applied to `res.samples` and `res.importance_weights()`). `DynamicNestedSampler` when the posterior matters more than the evidence. Evidences depend on prior volumes, so a Bayes factor is reported with both priors. 10. Summaries: median with 16th and 84th percentiles per parameter, a corner plot, posterior predictive curves over the data, and derived quantities computed from the samples rather than from point estimates. For simple propagation use `astropy.uncertainty.Distribution` or sample arithmetic; state correlations. ## Report - Data: source, filter or instrument, N points, time scale and reference, quality cuts. - Model, parameterization, prior table, sampler settings (walkers, steps, nlive, seed). - Diagnostics: tau per parameter, N / tau, acceptance fraction, R-hat, dlogz reached. - Parameter table with credible intervals, evidence if computed, and known systematics. ## Sources - astropy documentation (units, coordinates, time, FITS, WCS): https://docs.astropy.org/en/stable/ - astropy Lomb-Scargle periodograms: https://docs.astropy.org/en/stable/timeseries/lombscargle.html - emcee autocorrelation analysis and convergence: https://emcee.readthedocs.io/en/stable/tutorials/autocorr/ - dynesty documentation: https://dynesty.readthedocs.io/en/stable/ - ArviZ diagnostics (`rhat`, `ess`): https://python.arviz.org/en/stable/api/diagnostics.html - VanderPlas (2018), Understanding the Lomb-Scargle periodogram, ApJS 236, 16: https://doi.org/10.3847/1538-4365/aab766 - Eastman, Siverd and Gaudi (2010), Achieving better than 1 minute accuracy in the HJD and BJD, PASP 122, 935: https://doi.org/10.1086/655938 - Hogg and Foreman-Mackey (2018), Data analysis recipes, using MCMC, ApJS 236, 11: https://doi.org/10.3847/1538-4365/aab76e