SDSS fit: RA=146.4905, Dec=-0.3693

This notebook fetches an SDSS spectrum at fixed sky coordinates and fits with fit_pl=True, fit_fe=False, fit_bc=False.

Current default: fit_method="optax" runs staged SVI/Optax MAP optimization (continuum warm start, then full model). fit_method="optax+nuts" uses that staged MAP point to initialize NUTS.

fit_poly_edge_flex is no longer used; continuum flexibility is controlled by fit_poly=True and fit_poly_order.

import numpy as np
from astroquery.sdss import SDSS
from astropy.coordinates import SkyCoord
from astropy import units as u

import jaxqsofit
from jaxqsofit import QSOFit

coord = SkyCoord(146.4905, -0.3693, unit="deg")
xid = SDSS.query_region(coord, spectro=True, radius=5 * u.arcsec)
sp = SDSS.get_spectra(matches=xid[:1])[0]

tb = sp[1].data
lam = np.asarray(10 ** tb["loglam"], dtype=float)
flux = np.asarray(tb["flux"], dtype=float)
ivar = np.asarray(tb["ivar"], dtype=float)

err = np.full_like(flux, 1e-6)
m = np.isfinite(ivar) & (ivar > 0)
err[m] = 1.0 / np.sqrt(ivar[m])

z = float(sp[2].data["z"][0])
ra = float(coord.ra.deg)
dec = float(coord.dec.deg)
plateid = int(sp[0].header.get("plateid", 0))
mjd = int(sp[0].header.get("mjd", 0))
fiberid = int(sp[0].header.get("fiberid", 0))
sdss_filename = f"{plateid:04d}-{mjd}-{fiberid:04d}"

lam.size, z, sdss_filename

q = QSOFit(
    lam=lam,
    flux=flux,
    err=err,
    z=z,
    ra=ra,
    dec=dec,
    filename=sdss_filename,
    output_path=".",
)

q.fit(
    deredden=True,
    fit_method="optax+nuts",
    fit_lines=True,
    decompose_host=True,
    fit_pl=True,
    fit_fe=False,
    fit_bc=False,
    fit_poly=True,
    fit_poly_order=2,
    dsps_ssp_fn="../tempdata.h5",
    fsps_age_grid = (0.01, 0.03, 0.1, 0.3, 1.0, 3.0, 10.0),
    fsps_logzsol_grid = (-1.0, -0.5, 0.0, 0.2),
    optax_steps=600,
    optax_lr=1e-2,
    nuts_warmup=50,
    nuts_samples=50,
    nuts_chains=1,
    plot_fig=True,
    save_fig=False,
    save_result=False,
)

s = q.numpyro_samples

# Blue-side RBF width
m16, m50, m84 = np.nanpercentile(np.asarray(s["edge_rbf_sigma_blue"], dtype=float), [16, 50, 84])
print(f"edge_rbf_sigma_blue = {m50:.3f} (+{m84-m50:.3f}/-{m50-m16:.3f})")

# Blue-side RBF amplitudes (vector: one per basis function)
amps_blue = np.asarray(s["edge_rbf_amp_blue"], dtype=float)  # shape: [nsamp, n_blue_basis]
for j in range(amps_blue.shape[1]):
    m16, m50, m84 = np.nanpercentile(amps_blue[:, j], [16, 50, 84])
    print(f"edge_rbf_amp_blue[{j}] = {m50:.3f} (+{m84-m50:.3f}/-{m50-m16:.3f})")

s = q.numpyro_samples
pl_norm = np.asarray(s["PL_norm_eff"], dtype=float)
pl_slope = np.asarray(s["PL_slope"], dtype=float)
pivot_wave = float(getattr(q, "pivot_wave", 0.5 * (q.wave[0] + q.wave[-1])))

f_lambda_2500 = pl_norm * (2500.0 / pivot_wave) ** pl_slope

c_A_s = 2.99792458e18
f_nu = (f_lambda_2500 * 1e-17) * (2500.0**2) / c_A_s
m_2500_samples = -2.5 * np.log10(f_nu) - 48.60

m16, m50, m84 = np.nanpercentile(m_2500_samples, [16, 50, 84])
print(f"pivot_wave = {pivot_wave:.1f} A")
print(f"m_2500 = {m50:.3f} (+{m84-m50:.3f}/-{m50-m16:.3f})")

Estimate Host and BC fractions

In this case, the host fraction is an upper-limit because the PL has a tendency to soak up things even when no AGN is there.

frac_host_samp = 1.0 / (1.0 + np.exp(-np.asarray(s['log_frac_host'])))
r16, r50, r84 = np.nanpercentile(frac_host_samp, [16, 50, 84])
print(f"host fraction = {r50:.3f} (+{r84-r50:.3f}/-{r50-r16:.3f})")

i3000 = np.argmin(np.abs(np.asarray(q.wave) - 3000.0))

bc_draws = np.asarray(q.pred_out["f_bc_model"], dtype=float)[:, i3000]
pl_draws = np.asarray(q.pred_out["f_pl_model"], dtype=float)[:, i3000]

bc_over_pl_draws = bc_draws / pl_draws

r16, r50, r84 = np.nanpercentile(bc_over_pl_draws, [16, 50, 84])
print("lambda =", q.wave[i3000])
print(f"BC/PL at ~3000A = {r50:.3f} (+{r84-r50:.3f}/-{r50-r16:.3f})")