A scalable Bayesian framework for galaxy emission line detection and redshift estimation

A SCALABLE B A YESIAN FRAMEWORK FOR GALAXY EMISSION LINE DETECTION AND REDSHIFT ESTIMA TION B Y A L E X A N D E R K U H N 1 , a , B O N NA B E L L E Z A B E L L E 2 , d , S A R A A L G E R I 1 , b , G A L I N L . J O N E S 1 , c , A N D C L AU D I A S C A R L A T A 2 , e 1 School of Statistics, University of Minnesota, a kuhn0218@umn.edu ; b salgeri@umn.edu ; c galin@umn.edu 2 School of Physics and Astr onomy , University of Minnesota, d kusch037@umn.edu ; e mscarlat@umn.edu Estimating galaxy redshifts is crucial for constraining ke y ph ysical quan- tities like those in the equation of state of dark energy . Modern telescopes such as the James W ebb Space T elescope, the Euclid Space T elescope, and the N ASA Nancy Grace Roman Space T elescope are producing massi ve amounts of spectroscopic data that enable precise redshift estimation. How- ev er, a galaxy’ s redshift can be estimated only when emission lines are present in the observed spectrum, which is unknown a priori . A novel Bayesian ap- proach to estimating redshift and simultaneously testing for the presence of emission lines is dev eloped. Although modern spectroscopic surveys in volv e millions of spectra and giv e rise to highly multimodal posterior distrib utions, the proposed framew ork remains computationally efﬁcient, admitting a par- allelizable implementation suitable for large-scale inference. 1. Introduction. Despite signiﬁcant progress in understanding the lar ge-scale structure and e volution of the universe, major uncertainties remain about the nature of dark ener gy and dark matter . These components continue to elude direct detection and challenge our understanding of fundamental physics. Se veral ongoing and planned experiments aim to address these uncertainties. The Euclid space telescope (a European Space Agency mission with contributions from N ASA) and the future N ASA Nancy Grace Roman Space T elescope are tw o such missions, designed to map the geometry of the dark uni verse and inv estigate the cause of its accelerating e xpansion. Both will construct three-dimensional maps of the univ erse by measuring galaxy redshifts , which serve as proxies for distance. As the univ erse expands, light from more distant galaxies is stretched to longer wav elengths, and redshift quantiﬁes this stretching. Because the precise relationship between distance and redshift depends on the underlying cosmological model, these measurements can place strong constraints on the nature of dark energy ( Laureijs et al. , 2011 ; Euclid Collaboration , 2025 ). High-precision galaxy redshifts are obtained from spectr oscopic data , which measure light intensity across wa velengths to produce spectra . Traditional slit-based spectroscopy works by isolating a single galaxy’ s light with a narro w slit and dispersing only that selected beam. T o observe many galaxies simultaneously , Euclid and Roman instead use slitless spectr oscopy , which disperses light from all galaxies in a telescope’ s ﬁeld of vie w at once, producing spectra for numerous galaxies in a single exposure. This high-throughput approach enables large surve ys, b ut introduces additional challenges for determining galaxy redshifts. Spectra can ov erlap, and the observed light at a giv en wa velength may include contributions unrelated to any particular galaxy . Slitless spectroscopy also suffers from a lower signal-to-noise ratio and reduced spectral resolution compared to slit-based methods ( Baronchelli et al. , 2020 ). W ithin these spectra, the most informativ e features for redshift estimation are bright peaks produced by elements such as hydrogen and oxygen, called emission lines . These lines hav e K e ywor ds and phrases: Bayesian inference, Signal detection, Slitless spectroscopy, Astrostatistics. 1 2 kno wn r est-frame wa velengths – the wa velengths they would occur at if the galaxy were stationary relati ve to the observer . In observed spectra, emission lines appear at higher w av e- lengths than their rest-frame values due to the galaxy’ s redshift, b ut otherwise retain their characteristic shapes (see Section 2.2 for details). Importantly , not ev ery spectrum contains emission lines, and their presence and observed locations are not known a priori . While e ven a single line can provide information about a galaxy’ s redshift, detecting multiple lines greatly improves precision, especially in the presence of spectral overlap and in cases with lo wer signal-to-noise ratios. A central aim of this work is to provide principled uncertainty quantiﬁcation for redshift estimation that accounts for this variability in emission-line pres- ence and strength. Roman and Euclid provide slitless spectroscopic capabilities ov er lar ge areas of the sky (on the order of half a square de gree per pointing). The James W ebb Space T elescope (JWST) of fers a complementary approach, performing slitless spectroscopy with its Near Infrared Im- ager and Slitless Spectrograph (NIRISS) and Near Infrared Camera (NIRCam) instruments. JWST covers smaller ﬁelds of view (a few thousandths of a square degree) at lower spectral resolution, b ut reaches substantially greater depths than Euclid and Roman. Despite these dif ferences in surv ey area and depth, the fundamental challenges in measuring spectroscopic redshifts remain similar . In this manuscript, we present an analysis of data from JWST as a testbed; ho we ver , the methodology is designed with the lar ge surve y volumes expected from Euclid and Roman in mind. Accordingly , our approach prioritizes computational ef ﬁciency to enable the analysis of millions of spectra. 1.1. Main challeng es and existing work. W e develop a Bayesian approach to simultane- ously detect multiple emission lines and estimate redshift that explicitly accommodates the most common systematic errors induced by slitless spectroscopic data. T o the best of the authors’ kno wledge, no other existing approach has done so. Most redshift estimation pipelines currently rely on template ﬁtting , typically through cross-correlation maximization ( T onry and Davis , 1979 ; Zamora and Díaz , 2024 ) or χ 2 - minimization ( Bolton et al. , 2012 ). These methods match the observed spectrum to astro- physical templates constructed from simulations, high-quality observations, or both. While simple and f ast, template-ﬁtting methods do not pro vide a principled frame work for emission line detection or for quantifying redshift uncertainty . In an attempt to address the limitations of traditional template ﬁtting, Jamal et al. ( 2018 ) construct a Bayesian model based on a set of known templates and utilize Bayes factors to test for the presence of emission lines. This provides posterior redshift estimates and a princi- pled method for line detection when the a vailable set of templates adequately represents the observed spectra. Howe ver , a relev ant template may not exist for ev ery spectrum, introduc- ing unaccounted model uncertainty and increasing the risk of false line detections. Our work addresses these issues by using a template-free model. Other template-free approaches exist. Howe ver , they are limited to searching for a sin- gle emission line rather than multiple lines simultaneously . These include machine learning- based classiﬁcation methods ( Baronchelli et al. , 2020 , 2021 ) and hierarchical Bayesian mod- els ( Park, V an Dyk and Siemiginowska , 2008 ; Harrison, Lochner and Brown , 2017 ). These Bayesian approaches typically rely on Markov Chain Monte Carlo (MCMC) for posterior inference, which can be especially challenging in this setting. Even when searching for a single emission line, the posterior can be multimodal, as illustrated by Park, V an Dyk and Siemigino wska ( 2008 ) in X-ray spectroscopy . When considering multiple lines simultane- ously , the likelihood is non-smooth and more sev erely multimodal due to the changing set of observ able emission lines for different redshift v alues: a challenge not addressed in previous work. B A YESIAN EMISSION LINE DETECTION 3 0 1 2 3 4 0.00 0.05 0.10 0.15 0.20 Redshift Marginal posterior density Strong line No line F I G 1 . Discrete approximation to the posterior mar ginal density (grid r esolution of 0.001) for two simulated datasets, one with a str ong line pr esent, and the other with no lines present. Both ar e normalized so that the r esulting pr obability mass functions sum to one. Figure 1 shows a discrete approximation to the marginal posterior of redshift (described in Section 5 ) when a strong signal is present, and when no signal is present. In the no-signal case, the estimated posterior is diffuse, placing small probabilities on many redshift values. In the strong-signal case, the estimated posterior places higher mass on a much smaller set of redshift v alues; howe ver , the resulting posterior is still multimodal, with large separation be- tween modes. This illustrates a difﬁcult interaction between redshift and signal strength that makes jointly modeling the two especially computationally challenging. Thus, ev en sophisti- cated MCMC methods including tempering-based algorithms ( Geyer and Thompson , 1995 ) or those with deri vati ve-based proposals ( Roberts and Rosenthal , 1998 ), may mix slo wly and can tak e days to produce reliable posterior estimates. This issue is particularly rele vant to the current work, as we use Bayes factors for emission line detection, which typically require substantially longer MCMC runs than standard posterior estimation to achie ve stable esti- mates. T o address these challenges, we propose an approximate approach that ev aluates the posterior ov er a discrete grid of redshift v alues, lev eraging conditional conjugac y to explicitly decouple the redshift and signal intensities. This parallelizable procedure is fast and scales well for massiv e datasets, reducing computation time from hours or days (as with MCMC) to seconds. It is worth noting that taking a frequentist approach is dif ﬁcult in this setting, since the inference problem is irre gular . In particular , under the null model of no signal, redshift is not identiﬁable, and thus classical asymptotic results such as W ilk’ s theorem ( W ilks , 1938 ) do not apply . T o address this, early w ork by Davies ( 1977 , 1987 ) introduced a theoretical frame work for hypothesis testing based on a likelihood ratio process indexed by the non-identiﬁable pa- rameter . This line of work was later extended in the high-energy ph ysics literature to address the look-else wher e effect (LEE) – the inﬂation of false discovery probability that arises when searching for a signal over a region, rather than testing at a prespeciﬁed location ( Gross and V itells , 2010 ). Howe ver , existing frequentist methods that handle the LEE are not directly ap- plicable to the problem at hand as the required regularity conditions (see, e.g., Algeri and v an Dyk , 2021 ) are not met. For example, for each ﬁxed value of redshift, the number of poten- tially activ e emission lines (and thus the number of constraints) is changing. This causes the 4 1.0 1.2 1.4 1.6 1.8 2.0 2.2 −4 −2 0 2 4 6 8 10 Obser ved w avelength ( µ m) Flux (10 − 20 cgs / A ° ) F I G 2 . An observed galaxy spectrum fr om JWST . The data ar e binned corr esponding to the ﬂat r ectangular natur e of the plot. The gre y dotted lines repr esent + / − two standard deviations ar ound the observed ﬂux values using the r eported measur ement err ors fr om JWST . limiting process to change along with the value of redshift, which is a dif ﬁculty not accounted for in existing theoretical work. Furthermore, resampling-based approaches ( Hansen , 1996 ) that may avoid such issues are computationally infeasible in surve ys containing millions of spectra and a desired signiﬁcance lev el of 3 σ (i.e., le vel Pr( Z > 3) ≈ 0 . 00135 in the one- sided case, where Z ∼ N (0 , 1) ). An alternativ e approach frames the problem of detecting a signal with unkno wn location as a multiple testing problem, despite the search space being inherently continuous. This can be addressed by discretizing the space and applying, for instance, a Bonferroni correction. Ho wev er, this will be extremely conservati ve when the search space is large ( Algeri et al. , 2016 ) as is the case in the current application, where thousands of tests would need to be conducted. Finally , Bayer and Seljak ( 2020 ) attempt to address the LEE by using Bayes factors, but they employ improper priors, which lead to non-unique Bayes factors. The use of improper priors with Bayes factors requires care to av oid introducing arbitrary constant terms that can lead to unreliable inference ( Berger and Pericchi , 1996 ). This is the ﬁrst work to allow for simultaneous detection of multiple emission lines with- out reliance on templates, balancing the trade-of fs necessary to be both (i) ﬂexible enough to account for common systematic uncertainties in the model’ s background and noise compo- nents, and (ii) computationally fast so that the method is practical for massive data surve ys e ven in the presence of a challenging posterior distrib ution. 2. Data and physical properties. 2.1. Data description. W e analyze slitless spectroscopic data obtained with JWST’ s NIRISS and NIRCam instruments. Because slitless spectroscopy disperses light from many galaxies simultaneously , each detector pixel records photons that may originate from multiple galaxies. Data-reduction pipelines then separate these o verlapping contrib utions to produce a one-dimensional spectrum for each galaxy . Since each pixel accumulates a large number of B A YESIAN EMISSION LINE DETECTION 5 OII Hg OIII Ha SII [SIII] [SIII] P ab OII Hg OIII Ha SII Redshift: 0 Redshift: 2 0.5 1.0 1.5 2.0 Obser ved w avelength ( µ m) F I G 3 . A comparison of observed emission line locations at r edshift 0 (top) and r edshift 2 (bottom). The light blue r ectangles r epr esent the observable range of the JWST detector . The OII, OIII, and SII doublets ar e visualized as single lines for simplicity . photon counts (e.g., at least 500) and the observ ations undergo extensi ve calibration and pro- cessing in JWST’ s data-reduction pipeline, the resulting spectra are reported as real-v alued ﬂux measurements. These are continuous-v alued measurements of light intensity , accompa- nied by associated measurement uncertainties. Due to the detector resolution, the observed wa velength range, [1 . 03 µ m , 2 . 25 µ m] (mi- crons), is di vided into N bins of similar but unequal widths. Additionally , JWST spectral data contain two detector gaps as a direct by-product of the observational strategy . More specif- ically , in all of the JWST spectra under consideration, observ ations are missing for wav e- lengths [1 . 28 µ m , 1 . 33 µ m] and [1 . 67 µ m , 1 . 78 µ m] . Denote by X obs the collection of inter- v als of wa velengths that were observed. Then, for JWST , X obs = [1 . 03 µ m , 2 . 25 µ m] \ X miss , where X miss = [1 . 28 µ m , 1 . 33 µ m] ∪ [1 . 67 µ m , 1 . 78 µ m] . In our analysis, we consider 213 high-quality JWST spectra that experts ha ve ﬂagged as containing emission lines. These spectra hav e accompanying redshift estimates determined by astronomers using a procedure described in Huberty et al. ( 2026 ). Brieﬂy , an automated line-ﬁtting procedure is ﬁrst applied, assuming that the strongest emission line corresponds to the element H α . Experts then visually assess the resulting ﬁt and test alternativ e emission line identiﬁcations until a satisfactory ﬁt is obtained. An example of a high-quality JWST spectrum is shown in Figure 2 . In addition, we consider 2800 JWST spectra that were not ﬂagged as high quality and were not selected by experts as containing emission lines. These spectra were not accompanied by astronomer-provided redshift estimates in the dataset used in this study , so the analysis in Section 6 provides no vel redshift estimates for these galaxies. 2.2. Emission lines and redshift. W e are interested in the possible detection of 13 emis- sion lines from various elements, each with a kno wn lab-frame wa velength (in µ m): OII (0.37271, 0.37299), H γ (0.434), H β (0.486), OIII (0.496, 0.501), H α (0.656), NII (0.658), SII (0.672, 0.673), [SIII] (0.907, 0.953), and Pa β (1.282). Some elements produce multiple closely spaced lines, kno wn as doublets , corresponding to the pairs of wa velengths listed abov e. The intensities of lines within each doublet are related by known constant factors, 6 follo wing well-established physical relationships (see Section 3.2 ). Similarly , the hydrogen lines (H γ , H β , and H α ) hav e a kno wn ordering in their relati ve strengths, with the exact v alues depending on the amount of dust obscuration present in each galaxy (see Section 4.2 ). When a galaxy is observed, each lab-frame wa velength, x lab k , is shifted by the unkno wn redshift, ζ , so that the observed wavelength is gi ven by x lab k (1 + ζ ) . This deterministic scaling determines which emission lines fall within the observ able wa velength region, X obs , linking observed emission line location to galaxy redshift. Figure 3 illustrates how the emission lines move into and out of the observable range as a function of redshift. This highlights why simultaneous detection of multiple lines is nec- essary: some lines may be unobservable at certain redshifts, while others remain detectable, providing complementary information. 3. Model likelihood. Assume the ﬂux values, denoted Y i for i = 1 , ..., N , are normally distributed. This section describes how the likelihood can be speciﬁed to incorporate not only the e xpected structure of the ﬂux v alues b ut also the rules governing emission-line eligibility , physical constraints, and the cov ariance structure of the observations. 3.1. Mean structur e. Each observ ation corresponds to the central wav elength of the i -th bin, denoted x i . The conditional mean of each Y i is modeled as (1) E [ Y i | α, f b , η , ζ , δ, x i ] = α + f b ( x i ) + 13 X k =1 η k ( ζ ) s k ( x i , ζ , δ ) , where α is an intercept term, f b ( · ) is an unknown background function, and the sum repre- sents the signal contribution from all 13 emission lines of interest. The background function, f b ( · ) , is expected to vary smoothly ov er the wavelength range, although little additional structure is known a priori . The limited kno wledge about f b ( · ) arises from the data processing required to produce the spectra and from possible contami- nants (“glitches”) in the data. The background is modeled as f b ( x ) = J X j =1 β j b j ( x ) , where b j ( · ) denotes the j -th cubic B-spline basis function. W e ﬁx J = 15 for the analysis in Section 6 . Sensiti vity to this choice is examined in Section S.7.2 of the Supplementary Material and summarized in Section 7 . The signal component is a linear combination of emission line intensities, η k ( ζ ) , and kno wn line-shape functions, s k ( · , ζ , δ ) , which depend on wav elength, the redshift parameter , ζ , and the line-width parameter , δ . Each line-shape function is modeled as a Gaussian density centered at the redshifted lab-frame wav elength and ev aluated ov er the observed wa velength range x ∈ X obs : (2) s k ( x, ζ , δ ) = 1 √ 2 π δ 2 w 2 exp  − 1 2 δ 2 w 2 [ x − x lab k (1 + ζ )] 2  , k = 1 , ..., K , where w is the median pixel width. The line-width parameter , δ , scales the width of each emission line width in units of w . 3.2. Emission line eligibility and constraints. Let η k ( ζ ) denote the signal intensity cor- responding to the k -th emission line at redshift ζ , and deﬁne the column v ector of all 13 signal intensities as η ( ζ ) = ( η 1 ( ζ ) , . . . , η 13 ( ζ )) ′ B A YESIAN EMISSION LINE DETECTION 7 where the prime symbol indicates the transpose. A line is considered eligible at redshift ζ if its redshifted lab-frame wavelength falls in the observ able region. Hence, the set of eligible indices is E ( ζ ) = { k : x lab k (1 + ζ ) ∈ X obs } . For k ∈ E ( ζ ) , we set η k ( ζ ) = 0 . Certain pairs of emission line intensities satisfy deterministic physical relationships, as de- scribed in Section 2.2 . Let C denote the set of such ordered pairs ( k , ℓ ) with k < ℓ . Whenev er both lines in a constrained pair are eligible, the corresponding signal intensities satisfy η ℓ ( ζ ) = c kℓ η k ( ζ ) , ( k , ℓ ) ∈ C ; k , ℓ ∈ E ( ζ ) . Requiring both lines to be eligible in order to impose the constraint is necessary here to av oid the case where one line falls in the observ able region, b ut the other does not. In that case, we still need to estimate the eligible line freely . The speciﬁc emission line constraints are (3) η 2 ( ζ ) = η 1 ( ζ ) , η 6 ( ζ ) = 0 . 336 η 5 ( ζ ) , η 10 ( ζ ) = 0 . 738 η 9 ( ζ ) , η 12 ( ζ ) = 0 . 4 η 11 ( ζ ) . T o av oid estimating redundant parameters, deﬁne K ( ζ ) ⊆ E ( ζ ) as the set of eligible indices that remain after removing the second inde x of each constrained pair when both members are eligible (so the smaller index, k , is retained by con vention). Formally , (4) K ( ζ ) = E ( ζ ) \ { ℓ : ( k , ℓ ) ∈ C , and k , ℓ ∈ E ( ζ ) } . The reduced signal intensity vector can no w be deﬁned as the subvector η sub ( ζ ) = { η k ( ζ ) } k ∈K ( ζ ) , and we only need to estimate η sub ( ζ ) . T o express the mean structure in terms of η sub ( ζ ) , begin by deﬁning the full signal matrix at ( ζ , δ ) as the matrix with elements: (5) [ W ( ζ , δ )] ik := s k ( x i , ζ , δ ) , i = 1 , . . . , N , k = 1 , . . . , 13 . Columns corresponding to deterministically related lines are combined through linear trans- formations so that only one coef ﬁcient in each constrained pair is estimated (again, using the smaller index by con vention). For example, if both η 5 and η 6 are eligible, their respective columns in W ( ζ , δ ) are mer ged into a single column with entries s 5 ( x i , ζ , δ ) + 0 . 336 s 6 ( x i , ζ , δ ) , i = 1 , . . . , N . Columns corresponding to ineligible lines are then remov ed, yielding the reduced signal ma- trix, W sub ( ζ , δ ) , of dimension N × |K ( ζ ) | , where |K ( ζ ) | denotes the cardinality of K ( ζ ) . Note that while removing columns from W ( ζ , δ ) and entries from η ( ζ ) corresponding to ineligible lines is not strictly necessary , reducing the dimension to obtain W sub ( ζ , δ ) leads to straightforward deri vations for the posterior distrib ution in terms of generalized least-squares estimates (see Sections S.2 and S.3 of the Supplementary Material). 3.3. Covariance structur e. The cov ariance of Y = ( Y 1 , . . . , Y N ) ′ is modeled as (6) S ( σ 2 ) = σ 2 I N + C , where I N is the N × N identity matrix, σ 2 is an unkno wn scalar , and C is a kno wn diagonal matrix with entries C ii = c 2 i gi ven by the reported measurement errors. In practice, the matrix C is usually diagonal, but the methodology remains unchanged if it were not. The additiv e structure in equation ( 6 ) allows explicit incorporation of measurement error information from JWST through the matrix C , which is necessary to pre vent false detec- tions. Measurement errors can become large near the edges of the observable region, X obs , producing extreme ﬂux values that might otherwise be mistaken for emission lines. The ho- moskedastic noise term, σ 2 I N , is included to capture additional statistical v ariability . 8 3.4. Linear model formulation. W e can now express the model for the ﬂux as a linear model for ﬁxed redshift and line-width pairs, ( ζ , δ ) , as (7) Y = α 1 N + X β + W sub ( ζ , δ ) η sub ( ζ ) + ε, where 1 N is the length- N vector of ones, X is a matrix of elements X ij = b j ( x i ) for i = 1 , . . . , N and j = 1 , . . . , J , and β = ( β 1 , . . . , β J ) ′ . Additionally , assume ε = ( ε 1 , . . . , ε N ) ′ ∼ N (0 , S ( σ 2 )) . T o simplify posterior deriv ations, we use orthogonalized versions of X and W sub ( ζ , δ ) with respect to S ( σ 2 ) deﬁned in Section S.1 of the Supplementary Material. 4. Priors. W e consider two competing models for the observed spectrum. Under the null model, M 0 , no emission lines are present. Under the alternativ e model, M 1 , one or more emission lines are present, with unknown redshift ζ , line width δ , and signal intensities η sub ( ζ ) . The background parameters α and β , as well as the noise variance σ 2 , are common to both models. In this section, we specify prior distributions for all parameters appearing under each model. For the full parameter vector , θ = ( α, β , σ 2 , η ( ζ ) , ζ , δ ) ′ , we assume the following prior structure (8) p ( θ ) ∝ p ( α ) p ( β ) p ( σ 2 ) p ( η sub ( ζ ) | ζ ) p ( ζ ) p ( δ ) , where p ( · ) denotes a probability density function with respect to the appropriate dominating measure. Here, η sub ( ζ ) is deﬁned as in Section 3 and contains the free elements of η ( ζ ) at redshift ζ . The remaining elements in η ( ζ ) are deterministically speciﬁed by the redshift and η sub ( ζ ) . Therefore, they affect the prior only by a constant factor . Priors for the parameters contributing to the product on the right-hand side of equation ( 8 ) are speciﬁed in the following subsections, along with hyperparameter choices motiv ated by existing astrophysical theory and expert input. A sensitivity analysis is presented in Section S.7 of the Supplementary Material and summarized in Section 7 . 4.1. Priors common to both models. For the intercept term, we assume α ∼ N ( a 0 , b 2 0 ) with a 0 = 5 and b 0 = 10 , reﬂecting the prior belief that ﬂux values are typically nonnega- ti ve, while allo wing for negati ve reported ﬂux values arising from the JWST pre-processing pipeline. For the B-spline background coefﬁcients, we assume β j ∼ N ( a 1 , b 2 1 ) , with a 1 = 0 and b 1 = 30 , for j = 1 , . . . , 15 . The prior mean of zero reﬂects the initial best guess of a constant background, while the large prior variance and relati vely lar ge number of basis functions ( J = 15 ) allow ﬂexibility to capture contaminant and spectral o verlap effects. W e ﬁx the noise variance, σ 2 , at its maximum-likelihood estimate under M 1 , denoted ˆ σ 2 (see Section S.4 of the Supplementary Material), and use this ﬁxed v alue in both M 0 and M 1 for all subsequent inference. Although a conjugate in verse-Gamma prior on σ 2 would be computationally con venient, it would require inducing an e xplicit dependence of the prior v ariances for α, β , and η sub ( ζ ) on σ 2 , which is not physically meaningful in this setting. Non-conjugate alternativ es that integrate ov er σ 2 av oid this dependence, but are computa- tionally costly since σ 2 enters the cov ariance structure and would require repeated in versions of relati vely lar ge matrices. 4.2. Model speciﬁc priors. W e model the absence of a signal (the null model) via a point mass prior at η ( ζ ) ≡ 0 . Since no emission lines are present under this model, no prior for the redshift, ζ , is required. Under the alternati ve model, in which at least one line is present, we place a scaled Beta(3, 3) prior on redshift, ζ , o ver the interv al (0 , 4] , yielding a symmetric, parabolic density peaked B A YESIAN EMISSION LINE DETECTION 9 at ζ = 2 . This reﬂects a weak prior preference for the OII and H α lines abo ve, since the prior mean is set near their eligible range. The line-width parameter , δ , is given a uniform prior over [1 , 3] , reﬂecting the expected range of narrow emission line widths relativ e to the data resolution (recall that the Gaussian line variance is δ 2 w 2 , where w is the median pixel width). In practice, we discretize both ζ and δ to enable an efﬁcient computation strate gy described in Section 5 below . The signal intensities, conditional on ζ , are giv en the following product prior: p ( η sub ( ζ ) | ζ ) ∝ 1 { η sub ( ζ ) ∈A ( ζ ) } Y k ∈K ( ζ ) p ( η k ( ζ )) , where p ( η k ( ζ )) denotes a density with support on [0 , ∞ ) , and we deﬁne (9) A ( ζ ) :=  η sub ( ζ ) : η 4 ( ζ ) η 3 ( ζ ) ≥ 2 . 174 , η 7 ( ζ ) η 4 ( ζ ) ≥ 2 . 86 , η k ( ζ ) ≥ 0 , k ∈ K ( ζ )  . The set A ( ζ ) encodes the physical inequality constraints between the hydrogen line intensi- ties mentioned in Section 2.2 . Each p ( η k ( ζ )) is taken to be the density of a normal distribution with parameters a 2 ,k and b 2 ,k , leading to the multi variate truncated normal density: (10) p ( η sub ( ζ ) | ζ ) ∝ 1 { η sub ( ζ ) ∈A ( ζ ) } Y k ∈K ( ζ ) exp ( − 1 2 b 2 2 ,k ( η k ( ζ ) − a 2 ,k ) 2 ) , where the hyperparameter vectors (11) a 2 = ( a 2 , 1 , . . . , a 2 , 13 ) ′ , b 2 = ( b 2 , 1 , . . . , b 2 , 13 ) ′ collect the hyperparameters for all 13 emission lines. The multiv ariate truncated normal prior is con venient because it directly enforces linear inequality constraints among the hydrogen emission line intensities, while preserving some analytic tractability , thereby minimizing computational costs (see Section 5 ). Other prior families require constrained optimization or specialized sampling, which can substantially increase computation time. The speciﬁc v alues of a 2 and b 2 used for JWST are a 2 = (0 . 02 , 0 . 02 , 0 . 0092 , 0 . 02 , 0 . 2 , 0 . 596 , 0 . 2 , 0 . 02 , 0 . 02 , 0 . 0271 , 0 . 02 , 0 . 05 , 0 . 02) ′ , and b 2 = (0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 1 , 1 , 1 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01) ′ . For most emission lines, the prior mean is quite small, with inequality constraints in ( 9 ) imposed as equality constraints on the appropriate entries of a 2 . Additionally , when strong emission lines are present, they are more likely to be either OII or H α . This is reﬂected in the prior means, two orders of magnitude lar ger than the rest. Finally , for JWST , the model priors are p ( M 0 ) = 0 . 9 and p ( M 1 ) = 0 . 1 , reﬂecting prior belief that galaxies without observ able emission lines are much more common than those with observ able emission lines. 5. Model comparison and estimation. For signal detection, we compare the no-signal model, M 0 , with the signal model, M 1 , via their posterior odds, (12) p ( M 1 | y ) p ( M 0 | y ) = p ( y | ˆ σ 2 , M 1 ) p ( y | ˆ σ 2 , M 0 ) × p ( M 1 ) p ( M 0 ) , where p ( y | ˆ σ 2 , M h ) is the marginal density of the observed data, y , under model M h , h = 0 , 1 . The ﬁrst term in the product on the right-hand side of ( 12 ) is the Bayes factor , BF 10 , and the second term is the prior odds. The no-signal model, M 0 , includes the intercept and back- ground parameters, α and β , as well as the variance σ 2 . The signal model, M 1 , additionally includes the free signal intensities, η sub ( ζ ) , the redshift, ζ , and the line-width parameter , δ . 10 5.1. Mar ginal density under the null model. T o compute the marginal density under M 0 , recall that Y = α 1 N + X β + ε , where α ∼ N ( a 0 , b 2 0 ) , β ∼ N (0 , b 2 1 I ) , and ε ∼ N (0 , S ( σ 2 )) . Because Y | ˆ σ 2 , M 0 is a linear combination of normally distributed random variables, it fol- lo ws that Y | ˆ σ 2 , M 0 ∼ N ( a 0 1 N , S ( ˆ σ 2 ) + b 2 0 1 N 1 ′ N + b 2 1 X X ′ ) , where S ( ˆ σ 2 ) is deﬁned as in equation ( 6 ) with ˆ σ 2 in place of σ 2 . 5.2. Mar ginal density under the alternative model. Under M 1 , in which at least one emission line is present, the goal is to compute p ( y | ˆ σ 2 , M 1 ) based on the prior on η sub ( ζ ) | ζ in equation ( 10 ). Directly integrating over η sub ( ζ ) , α, β , ζ , and δ is intractable, so we con- sider a ﬁne grid of ( ζ , δ ) pairs and proceed as follows. Let ζ 1 , . . . , ζ R denote a predeﬁned, equally spaced grid ov er (0 , 4] , chosen to match the desired redshift resolution. Similarly , let δ 1 , . . . , δ D denote a predeﬁned grid ov er [1 , 3] . All calculations for M 1 are performed ov er the R × D grid of ( ζ r , δ ℓ ) pairs. 5.2.1. Inte gration over η sub ( ζ ) , α , and β . Consider a ﬁxed ( ζ r , δ ℓ ) pair . Rather than working with the truncated multiv ariate normal prior on η sub ( ζ r ) | ζ r in equation ( 10 ) di- rectly , we introduce an encompassing model in which the truncated prior is replaced by its unconstrained multi variate normal counterpart. Speciﬁcally , we deﬁne (13) η ∗ sub ( ζ r ) | ζ r ∼ N ( a 2 ( ζ r ) , B 2 ( ζ r )) , where η ∗ sub ( ζ r ) denotes the unconstrained v ersion of η sub ( ζ r ) , and (14) a 2 ( ζ r ) = { a 2 ,k } k ∈K ( ζ r ) , B 2 ( ζ r ) = diag ( { b 2 2 ,k } k ∈K ( ζ r ) ) denote the corresponding subv ector and diagonal submatrix of hyperparameters deﬁned in equation ( 11 ) for lines with indices in K ( ζ r ) . Let p ∗ ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) denote the density obtained by inte grating out α, β , and η ∗ sub ( ζ r ) . More explicitly , this density is given by (15) p ∗ ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) = Z p ( y | ˆ σ 2 , α , β , η ∗ sub ( ζ r ) , ζ r , δ ℓ ) × p ( α ) p ( β ) p ( η ∗ sub ( ζ r ) | ζ r ) dα dβ dη ∗ sub ( ζ r ) , where p ( y | ˆ σ 2 , α , β , η ∗ sub ( ζ r ) , ζ r , δ ℓ ) is the conditional density of y deﬁned by equation ( 7 ) (see Section S.1 of the Supplementary Material for an explicit expression). As in the null model, the conjugate Gaussian priors imply that p ∗ ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) is the density of a multi variate normal with mean a 0 1 N + W sub ( ζ r , δ ℓ ) a 2 ( ζ r ) and co variance matrix b 2 0 1 N 1 ′ N + b 2 1 X X ′ + W sub ( ζ r , δ ℓ ) B 2 ( ζ r ) W sub ( ζ r , δ ℓ ) ′ + S ( ˆ σ 2 ) . A routine calculation (see Section S.2.1 of the Supplementary Material) leads to (16) p ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) = p ∗ ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) Pr( η ∗ sub ( ζ r ) ∈ A ( ζ r ) | ˆ σ 2 , ζ r , δ ℓ , y , M 1 ) Pr( η ∗ sub ( ζ r ) ∈ A ( ζ r ) | ζ r , M 1 ) , where Pr( η ∗ sub ( ζ r ) ∈ A ( ζ r ) | ζ r , M 1 ) and Pr( η ∗ sub ( ζ r ) ∈ A ( ζ r ) | ˆ σ 2 , ζ r , δ ℓ , y , M 1 ) are the prior and posterior probabilities that η ∗ sub ( ζ r ) satisﬁes the constraints in equation ( 9 ). Details for computing these probabilities are provided in Section S.2 of the Supplementary Material; in particular, each probability reduces to a nonnegati vity constraint on a known multiv ariate normal distribution, which can be computed ef ﬁciently . B A YESIAN EMISSION LINE DETECTION 11 5.2.2. Inte gration over δ and ζ . Since δ has been discretized, integrating out δ for each ﬁxed ζ r is approximated by a Riemann sum ov er its support: p ( y | ˆ σ 2 , ζ r , M 1 ) = D X ℓ =1 p ( y | ˆ σ 2 , ζ r , δ ℓ , M 1 ) p ( δ ℓ )∆ δ , where ∆ δ denotes the grid spacing for δ (with ∆ δ = 0 . 5 in practice). The marginal density under the alternati ve model is then approximated by (17) p ( y | ˆ σ 2 , M 1 ) ≈ R X r =1 p ( y | ˆ σ 2 , ζ r , M 1 ) p ( ζ r )∆ ζ , where ∆ ζ denotes the grid spacing for ζ . For spectra observed by JWST , we use ∆ ζ = 0 . 001 . For the chosen priors, equation ( 16 ) must be ev aluated for all 6 × 4000 = 24 , 000 ( ζ r , δ ℓ ) pairs, but the computations can be easily parallelized. 5.3. Redshift estimation under the alternative model. The posterior mar ginal distribution ζ is approximated by (18) p ( ζ r | y , ˆ σ 2 , M 1 ) ≈ p ( y | ˆ σ 2 , ζ r , M 1 ) p ( ζ r )∆ ζ P R r =1 p ( y | ˆ σ 2 , ζ r , M 1 ) p ( ζ r )∆ ζ , for r = 1 , . . . , R . The discrete approximation in equation ( 18 ) is used directly for model selection and redshift estimation, since the grid resolution is suf ﬁciently ﬁne for inference. As a point estimate of the redshift, we use the approximate mar ginal maximum a posteriori (MAP) estimator , (19) ˆ ζ = arg max r p ( ζ r | y , ˆ σ 2 , M 1 ) . Point estimates for the remaining parameters are taken as their posterior means, conditional on ζ = ˆ ζ and σ 2 = ˆ σ 2 . Deri vations for these estimates are provided in Section S.3 of the Supplementary Material. Uncertainty in redshift estimation is summarized using approximate mar ginal highest pos- terior density (HPD) sets for ζ . Because the posterior, p ( ζ | ˆ σ 2 , y , M 1 ) , is ev aluated on a discrete grid and is often multimodal, the HPD region typically consists of multiple disjoint interv als rather than a single contiguous range. T o construct a (1 − α ) HPD set, we order the grid points { ζ r } R r =1 by decreasing posterior probability and include the most probable val- ues until the cumulativ e posterior mass reaches at least 1 − α . Adjacent included grid points are then grouped into contiguous components, each corresponding to a local mode of the posterior distrib ution. The union of these components forms the approximate (1 − α ) HPD set. 6. Results on JWST data. T o assess the performance of the proposed procedure on JWST data, we begin with three individual spectra that illustrate key components of the analysis. W e then present aggreg ate results for the 213 high-quality (labeled) spectra and compare our emission-line detection decisions and redshift estimates with those provided by JWST astronomers. Finally , we analyze 2800 lo wer-quality (unlabeled) JWST spectra that hav e not been assigned e xpert classiﬁcations, o wing to both the dif ﬁculty of characterizing these spectra and the size of the dataset. A subset of the JWST spectra (galaxies 00004, 00033, and 02263), along with associated astronomer-pro vided redshift annotations and R code required to reproduce these results, is provided in the Data and R Code Supplement. The full dataset is subject to JWST data-sharing restrictions but is a vailable upon request. 12 0 1 2 3 4 0.00 0.05 0.10 0.15 Redshift Marginal posterior density F I G 4 . Discr ete approximation to the posterior marginal density for spectrum 00004. The dashed orange line indicates the r edshift estimate pr ovided by astr onomers. Model ﬁt assessment and a sensiti vity analysis of prior hyperparameter choices, including an in vestigation of false positi ve rates, are presented in Sections S.6 and S.7 of the Supple- mentary Material. Emission line detection is based on the astronomers’ prior model proba- bilities, p ( M 0 ) = 0 . 9 and p ( M 1 ) = 0 . 1 , giv en in Section 4 . W e select M 1 if the log posterior odds are greater than zero, and M 0 otherwise. 6.1. Analyzing individual spectra. W e analyze spectra from three illustrative g alaxies in the high-quality (labeled) JWST sample, corresponding to galaxy IDs 00004, 00033, and 02263 in the Data and R Code Supplement. Spectrum 00004 has a relati vely constant back- ground and a single prominent Gaussian line in the observed range. Emission line detection is therefore straightforward, but redshift estimation is more difﬁcult, as the marginal poste- rior of redshift is multimodal. Spectrum 00033 has a more complex background. Ho we ver , multiple emission lines are detected across the observed range, which more tightly constrain the redshift. In this case, the marginal posterior is unimodal. Finally , spectrum 02263 has a prominent, non-constant background and relatively weak emission lines, making emission lines difﬁcult to identify . The e vidence for detection is weak (although emission lines are detected), and the redshift is only weakly constrained by the observed data. T ogether , these three cases demonstrate how background structure and the number of detectable emission lines inﬂuence both the e vidence for emission line detection and the shape of the redshift posterior . While these examples highlight recurring features in the data, the full sample ex- hibits a broader range of background shapes and emission line conﬁgurations. 6.1.1. Galaxy 00004. Consider the JWST spectrum shown in Figure 2 . The log poste- rior odds for M 1 versus M 0 is 26.47, corresponding to a posterior probability of M 1 being ef fectiv ely 1. For reference, the corresponding log Bayes factor is 28.67. These v alues pro- vide e xtremely strong e vidence for the presence of at least one emission line, consistent with astronomers’ high conﬁdence that this spectrum contains emission lines. Based on the strong e vidence for detection, we proceed to estimate the galaxy’ s red- shift. Figure 4 displays the discrete approximation to the posterior marginal distribution B A YESIAN EMISSION LINE DETECTION 13 −5 0 5 10 1.00 1.25 1.50 1.75 2.00 2.25 Obser ved w avelength ( µ m) Flux F I G 5 . A posterior pr edictive plot assessing model ﬁt for spectrum 00004 in Figur e 2 . The black dots show the observed ﬂux values per wavelength. Her e, 5000 replicate datasets ar e drawn fr om the posterior pr edictive distribution. The pointwise median per wavelength is shown in dark blue, and the pointwise 2.5% and 97.5% pointwise quantiles as the lower and upper bounds of the light blue r e gion, respectively . for redshift. The posterior is multimodal, with four modes (the fourth being quite small). The dominant mode is centered at 3.837, the marginal MAP for redshift. The correspond- ing marginal 99.865% HPD set (the 3 σ set, reﬂecting astrophysical con vention) is the union [0 . 890 , 0 . 894] ∪ [1 . 678 , 1 . 684] ∪ [1 . 744 , 1 . 749] ∪ [3 . 831 , 3 . 844] . This set contains the redshift v alue 1.74734 reported by astronomers, although the posterior mar ginal MAP differs substan- tially from their point estimate. Importantly , the astronomers’ redshift is pro vided without an accompanying measure of uncertainty and should not be regarded as the ground truth. In this case, the overlap between the HPD set and the astronomers’ estimate offers partial agreement, while the multimodality of the posterior highlights substantial uncertainty . Next, we report the set of non-zero emission lines at the mar ginal MAP redshift, and com- pare it to the emission lines implied by the astronomers’ point estimate. At redshift 3.837, the activ e lines include both OII lines and H γ . In this case, the large spike in Figure 2 corre- sponds to the OII doublet. The posterior means of the emission line intensities (conditional on ˆ σ 2 and ˆ ζ ) are 0.044, 0.044, and 0.007, respecti vely . In contrast, the astronomers’ reported redshift of 1.747 corresponds to emission lines of H γ , H β , the OIII doublet, H α , NII, and the SII doublet, with H α (and the nearby , weaker NII line) as the prominent line feature. This reﬂects the common modeling assumption that the strongest observed line corresponds to H α or OIII, without fully accounting for the possibility of other elements when emission lines are weak er . The prior incorporates this belief, b ut does not enforce it deterministically as in pre vious ﬁtting methods. Point estimates for the nuisance parameters α, β , σ 2 , and δ are computed as described in Section S.3 of the Supplementary Material, but are omitted here for bre vity . This spectrum illustrates that ev en when there is strong e vidence for the presence of at least one emission line, substantial uncertainty may remain in the posterior distribution for redshift. The full analysis required approximately 15 seconds on a 20-core processor . Next, we examine a posterior predicti ve plot to assess overall model ﬁt. The results are sho wn in Figure 5 . The black dots represent the observed ﬂux values. The light blue band corresponds to the pointwise 2.5% and 97.5% quantiles of 5000 draws from the posterior 14 10 20 30 1.00 1.25 1.50 1.75 2.00 2.25 Obser ved w avelength ( µ m) Flux F I G 6 . A posterior pr edictive plot assessing model ﬁt for spectrum 00033. The black dots show the observed ﬂux values per wavelength. Her e, 5000 r eplicate datasets fr om the posterior predictive distribution are drawn, with the pointwise median per wavelength shown in dark blue, and the pointwise 2.5% and 97.5% pointwise quantiles as the lower and upper bounds of the light blue r e gion, respectively . predicti ve distribution, and the dark blue line depicts the pointwise median of these draws. A full description of how posterior predictiv e samples are generated is provided in Section S.5 of the Supplementary Material. V isually , the model provides a good ﬁt to the data. The posterior predicti ve bands cap- ture most of the observed ﬂux values without being overly wide. In particular , the model accurately reproduces the sharp peak associated with the most prominent emission line. On the right-hand side of the plot, the predicti ve bands widen substantially near the boundary . This widening reﬂects the measurement errors provided by JWST astronomers, rather than instability in the ﬁtting procedure. 6.1.2. Galaxy 00033. W e no w analyze spectrum 00033 from the labeled JWST data. Figure 6 presents the corresponding posterior predicti ve plot. The log posterior odds in f av or of M 1 are 120.2749, indicating extremely strong evidence for the presence of at least one emission line. In contrast to spectrum 00004, the marginal 99.865% HPD set for redshift is the single interv al [1 . 867 , 1 . 871] , reﬂecting a unimodal marginal posterior , and the marginal MAP estimate is 1 . 869 . This is essentially identical to the astronomers’ reported redshift of 1 . 86884 . At this redshift, the detected emission lines are the OII doublet, H γ , H β , the OIII doublet, H α , NII, and the SII doublet. The corresponding posterior mean signal intensities (conditional on ˆ σ 2 and ˆ ζ ) are 0.0204, 0.0204, 0.0025, 0.0126, 0.0377, 0.1123, 0.0523, 0.0126, 0.0069, and 0.0093, respectiv ely , with the most prominent feature corresponding to the OIII doublet. The posterior predictiv e ﬁt is somewhat poorer than for spectrum 00004, as the predictiv e bands do not capture background and noise ﬂuctuations as closely . Ne vertheless, the pres- ence of multiple emission lines tightly constrains the redshift. In this case, strong detection e vidence coincides with precise redshift estimation. 6.1.3. Galaxy 02263. In contrast to the previously analyzed spectra, spectrum 02263 provides comparativ ely weak e vidence for emission-line detection due to a prominent as- trophysical background and faint emission lines. The posterior predictiv e plot is shown in B A YESIAN EMISSION LINE DETECTION 15 −10 0 10 20 1.00 1.25 1.50 1.75 2.00 2.25 Obser ved w avelength ( µ m) Flux F I G 7 . A posterior pr edictive plot assessing model ﬁt for spectrum 02263. The black dots show the observed ﬂux values per wavelength. Her e, 5000 r eplicate datasets fr om the posterior predictive distribution are drawn, with the pointwise median per wavelength shown in dark blue, and the pointwise 2.5% and 97.5% pointwise quantiles as the lower and upper bounds of the light blue r e gion, respectively . Log posterior odds Frequency 0 2 4 6 8 10 −10 10 30 50 70 90 110 140 170 200 230 260 290 F I G 8 . Histogr am of log posterior odds for 213 labeled spectra. Some values ar e cutoff for being too large . The dashed blue line at zer o repr esents the cutoff thr eshold for emission line detection. Figure 7 . For this spectrum, the log posterior odds in fav or of M 1 is 0.673, corresponding to a posterior probability of M 1 of 0.66. Thus, while at least one emission line is detected under the decision rule, the e vidence is substantially weaker than in the previous two cases. The marginal MAP for redshift is 1 . 286 , which is again quite similar to the astronomers’ provided redshift of 1 . 2854 . The corresponding 99.865% marginal HPD set for redshift the union of 13 disjoint intervals (omitted for brevity), reﬂecting considerable uncertainty in the redshift estimate. 16 213 astronomer redshifts Frequency 0 1 2 3 4 0 15 30 203 marginal MAP redshifts Redshift Frequency 0 1 2 3 4 0 15 30 0 1 2 3 4 0 1 2 3 4 Astronomer redshift Marginal MAP redshift F I G 9 . Left panel: Histograms of r edshift estimates for labeled spectr a fr om astr onomers (top) and marginal MAP estimates (bottom). Counts differ because our selection procedur e excludes ten spectra identiﬁed by astr onomers based on the log posterior odds. Right panel: Scatterplot comparing astr onomer and marginal MAP r edshift estimates for the 203 spectra selected by the log posterior odds. In terms of model ﬁt, the majority of observed ﬂux v alues fall within the pointwise pos- terior predictive intervals, except for a few outliers at lower wav elengths. The overall mean structure is reasonably well captured, and multiple weak emission features are visible in the posterior predicti ve median curv e. T aken together , these three examples demonstrate that redshift uncertainty is dri ven not only by the presence of emission lines but also by their number, relativ e strength, and inter- action with the background structure. A single prominent feature may yield decisi ve evidence for detection while leaving redshift weakly identiﬁed. In contrast, multiple strong lines can tightly constrain redshift. When the signal strength is weak relative to the background and noise lev els, both detection e vidence and redshift precision deteriorate. These regimes recur throughout the JWST dataset and are reﬂected more broadly in the aggregate analysis that follo ws. 6.2. Aggr e gate r esults for 213 labeled spectra. W e now consider the full set of 213 JWST spectra for which astronomers have identiﬁed emission lines and pro vided redshifts. Although these expert-provided annotations are not ground truth, the y provide a natural point of comparison for our emission line detection decisions and redshift estimates. Figure 8 displays a histogram of log posterior odds across the 213 spectra. The blue vertical line marks zero, the decision threshold for emission line detection. The majority of spectra lie to the right of this line, with large positi ve log-posterior odds in fa vor of M 1 . T en galaxy spectra are not selected as containing emission lines under the posterior odds decision rule. All 213 of these spectra were ﬂagged by astronomers as likely to contain emis- sion lines, so the high overall agreement is reassuring. Ne vertheless, our procedure is more conserv ativ e than the e xpert labeling in this subset. This difference reﬂects, in part, the use of prior model probabilities for a generic set of galaxy spectra rather than for a pre-selected set lik ely to contain emission lines. When using the log Bayes f actor alone (corresponding to equal prior model probabilities), only ﬁ ve spectra are not selected. B A YESIAN EMISSION LINE DETECTION 17 Labeled Unlabeled 0.001 0.005 0.050 0.500 T otal HPD width (log scale) Labeled Unlabeled 0 5 10 15 20 25 30 Number of inter vals in HPD set F I G 1 0 . Left panel: Boxplots of the total width of 99.865% HPD sets for spectra selected by our procedur e as containing at least one emission line for the labeled (left) and unlabeled (right) JWST datasets. HPD sets ar e generally much narr ower for the labeled data. Right panel: Boxplots of the number of disjoint intervals in each 99.865% HPD set. Labeled spectr a typically have 1–2 intervals, wher eas unlabeled spectr a mor e often have 3–6. Log posterior odds Frequency 0 100 200 300 400 500 600 −10 0 10 20 30 40 50 60 70 80 F I G 1 1 . Histogram of log posterior odds for 2715 unlabeled spectra. The dashed blue line at zero r epresents the cutoff decision for line detection. Note the r ange on the x-axis is smaller compar ed to Figur e 8 . Examination of the spectra not selected by our procedure indicates that they typically exhibit low signal-to-noise ratios or prominent background features relati ve to the emission line strengths (similar to spectrum 02263 abov e). In these cases, the redshift posterior is weakly constrained, so ev en a more liberal detection rule would not necessarily yield precise or informati ve redshift estimates. 18 Redshift Frequency 0 1 2 3 4 0 50 100 150 200 F I G 1 2 . Histo gram of r edshift mar ginal MAP estimates for spectra in the unlabeled set that ar e selected accor ding to posterior odds. 6.3. Aggr e gate r esults for 2800 unlabeled spectra. W e estimate the redshift only in spec- tra for which emission lines are detected. The left panel of Figure 9 compares histograms of the marginal MAP redshift estimates from our procedure with the redshifts reported by as- tronomers. W e emphasize again that neither set of values should be regarded as ground truth. Overall, the two distributions are similar . The most noticeable difference occurs in the low- redshift range: the astronomer-reported redshifts include no galaxies with redshifts between 0 and 0.5, whereas our procedure identiﬁes eight galaxies in this range. The right panel of Figure 9 shows a corresponding scatterplot comparing redshift estimates for the 203 spectra selected by our procedure based on the log posterior odds. There is strong agreement be- tween the two sets of estimates for most galaxies; howe ver , 22 galaxies lie of f the diagonal, indicating discrepancies in their ﬁtted redshifts. The left panel of Figure 10 summarizes the total widths of the mar ginal 99.865% HPD sets for redshift for spectra selected as containing emission lines by our procedure. The distribu- tion of total widths is highly skewed (note the log scale), with most labeled spectra having total widths between 0.005 and 0.02. For reference, recall that the redshift resolution used in the ﬁtting procedure is 0.001. T o better characterize the source of this uncertainty , the right panel of Figure 10 shows the number of disjoint interv als within each 99.865% HPD set. For the labeled data, most spectra hav e HPD sets consisting of a single interval, or occasionally two, indicating that the posterior is typically unimodal or bimodal. A small number of spectra exhibit many disjoint intervals, corresponding to more weakly constrained redshift estimates. Thus, for the labeled sample, the width of the 99.865% HPD sets generally reﬂects dispersion around a dominant redshift mode rather than the presence of se veral well-separated, competing redshift solutions. W e no w present aggreg ate results for the 2800 unlabeled JWST spectra. Due to the size of the dataset and the generally lower quality of these spectra, they hav e not been screened for emission line presence or assigned redshift estimates by astronomers. Within this set, some spectra exhibit data quality issues that render inference unreliable. These issues arise from various errors in the detection process, resulting in large numbers of e xtreme or missing v alues. The screening procedure used to select spectra of reasonable quality is described in B A YESIAN EMISSION LINE DETECTION 19 171 spectra from high−quality set Frequency 0 1 2 3 4 0 10 20 334 spectra from lo w−quality set Redshift Frequency 0 1 2 3 4 0 20 40 F I G 1 3 . Comparison of mar ginal MAP redshift distributions for labeled (top) and unlabeled (bottom) spectra after r estricting to spectra with str ong e vidence for emission line pr esence ( P ( M 1 | y ) > 0 . 99865 ), and 99.865% HPD sets consisting of one or two disjoint intervals. The dominant modes occur in dif ferent redshift r e gions acr oss the two samples. Section S.5 of the Supplementary Material and remov es 85 spectra from the analysis. W e therefore proceed using the remaining 2715 spectra. Figure 11 displays a histogram of log posterior odds for such spectra. The distribution is strongly right-ske wed, with most values belo w zero (more extreme v alues are truncated in the ﬁgure). The lo wer quartile, median, and upper quartile are -4.3, -2.6, and 4.1, respectiv ely , with a maximum v alue of 570. Under the posterior odds decision rule, emission lines are detected in 937 spectra, corresponding to approximately 35% of the unlabeled sample. Figure 12 displays a histogram of mar ginal MAP redshift estimates for the 937 spectra identiﬁed as containing at least one emission line. The ov erall shape of the distribution is broadly similar to that observ ed for the labeled sample. Howe ver , the modes near redshift 1.6 and 2.5 are more pronounced here than in the high-quality set. It is natural to ask whether spectra with particularly strong evidence for emission line pres- ence, or those with more tightly constrained redshift estimates (e.g., HPD sets consisting of a small number of disjoint interv als), produce a distribution more comparable to that of the high-quality spectra. T o examine this, we restrict attention to spectra with posterior proba- bility of M 1 at least 0.99865 (the 3 σ probability , corresponding to a log posterior odds of at least 6.6), and whose 99.865% HPD set for redshift consists of either one or two disjoint interv als. W e allo w up to two intervals in order to retain a reasonable number of observations from the high-quality sample, whose size is limited. The ﬁltered histograms for both the labeled and unlabeled spectra are sho wn in Figure 13 . The resulting distributions differ substantially . In particular , the dominant modes occur in distinct redshift regions across the two samples. Notably , the boundary redshift values ob- served in the high-quality set disappear entirely when restricting attention to spectra with strong e vidence for emission line presence. W e next examine posterior uncertainty in redshift estimation. Figure 10 summarizes pos- terior uncertainty for redshift in the unlabeled spectra (and labeled spectra as discussed in 20 the previous section). The total widths of the 99.865% HPD sets are generally larger than those observ ed in the labeled sample, and the number of disjoint interv als is correspondingly higher , indicating greater posterior uncertainty in redshift estimation. Nevertheless, 159 un- labeled spectra yield HPD sets consisting of a single interv al, suggesting well-constrained redshift estimates for a subset of galaxies that is nearly as large as the expert-labeled dataset. Moreov er , the posterior odds provide a natural ranking of spectra by the strength of evi- dence for emission-line detection. This ranking enables targeted follow-up analysis: spectra with large posterior odds can be prioritized for more sophisticated modeling where compu- tational cost is less restrictive, or post-hoc expert in vestigations. Processing all 2800 spectra required under 12 hours using modest computational resources (a 20-core processor and 32 GB of RAM). Because the ﬁtting procedure admits parallel computation within each spec- trum, the approach scales naturally to much larger datasets. 7. Final Remarks. W e ha ve presented a modeling framew ork and analysis of approx- imately 3,000 JWST spectra, providing a structured procedure for simultaneously detecting multiple emission lines and quantifying the uncertainty in redshift estimation. The approach explicitly accounts for systematic uncertainties arising from background modeling, excess v ariance beyond the reported measurement error , and emission line widths. Our analysis sug- gests that redshift uncertainty may be substantially larger than is typically reported, even for some high-quality spectra with expert labels. The results exhibit sensiti vity to the choice of some prior hyperparameters. In particular , the number of background basis functions J should be chosen sufﬁciently large to allow the null model adequate ﬂexibility (see Section S.7.2 of the Supplementary Material). The value J = 15 used in Section 6 appears appropriate for this purpose. Additionally , the log poste- rior odds and the distribution of mar ginal MAP redshift estimates are sensitiv e to the prior speciﬁcation for the signal intensities (Section S.7.3 ). Priors that incorporate astrophysical kno wledge of relati ve emission line strengths lead to more stable redshift inference, whereas dif fuse intensity priors can substantially alter the marginal MAP redshift distrib ution. This sensiti vity is ampliﬁed by using the marginal MAP as a point estimate, which depends on the relati ve heights of modes in the mar ginal posterior . In contrast, the results sho w relati vely little sensiti vity to the hyperparameters gov erning the redshift prior (Section S.7.4 ). W e recommend reporting HPD sets alongside marginal MAP estimates in downstream analyses whenev er possible. By capturing posterior multimodality , HPD sets are generally less sensiti ve to prior hyperparameter choices than the mar ginal MAP estimate alone. T o achiev e computational feasibility at the scale required by modern spectroscopic surv eys such as JWST , Euclid, and Roman, sev eral modeling compromises were necessary . These include distrib utional assumptions in the likelihood (e.g., normally distributed errors), the use of truncated multiv ariate normal priors for signal intensities, and a plug-in estimate for σ 2 in the co variance structure. These simpliﬁcations enable lar ge-scale analysis, b ut limit full Bayesian uncertainty propagation. Nev ertheless, the proposed methodology appears to perform well in our data analysis in Section 6 and in the numerical experiments (Sections S.6 and S.7 of the Supplementary Material). For analyses focused on speciﬁc galaxies, or reﬁned subsets where highly precise red- shift estimates are required, a fully Bayesian treatment would be desirable. Ev en for a single spectrum, howe ver , the posterior distribution can be highly multimodal. As illustrated in our analysis, a prominent emission feature may induce multiple competing redshift modes. Con- sequently , advanced MCMC methods such as parallel tempering and related tempered ap- proaches ( Ge yer and Thompson , 1995 ; Syed et al. , 2022 ) would likely be required. Dev elop- ing such algorithms for models with v arying numbers of acti ve emission lines and associated constraints remains an important methodological challenge. In addition, complementary fre- quentist procedures for emission line detection and redshift estimation would provide useful alternati ve inferential tools for lar ge spectroscopic surve ys. B A YESIAN EMISSION LINE DETECTION 21 Funding. AK, SA, and GJ were partially supported by NSF grant DMS-2152746. SA was partially supported by the W arwick Mid-Career A ward, CLA, Uni versity of Minnesota. AK was partially supported by the NSF Grant NR T -1922512. CS knowledges the support of N ASA ROSES Grant 12-EUCLID11-0004. SUPPLEMENT AR Y MA TERIAL Supplement to “A scalable Bayesian framework f or Galaxy emission line detection and redshift estimation” Additional deri vations, descriptions of optimization and sampling schemes, model ﬁt assess- ments, and a sensiti vity analysis. Data and R Code A subset of the full JWST spectral data, corresponding astronomer redshift annotations, and R code required to reproduce all results and ﬁgures reported in the paper . Materials also av ailable at https://github .com/alexkuhn0115/Bayesian- emission- line- detection- public . REFERENCES A L G E R I , S . and V A N D Y K , D . (2021). T esting one hypothesis multiple times. Statistica Sinica 31 959–979. https://doi.org/10.48550/arXi v .1701.06820 A L G E R I , S ., V A N D Y K , D . , C O N R A D , J . and A N D E R S O N , B . (2016). On methods for correcting for the look- elsewhere ef fect in searches for new physics. Journal of Instrumentation 11 P12010. https://doi.org/10.1088/ 1748- 0221/11/12/P12010 B A RO N C H E L L I , I . , S C A R L AT A , C . M ., R O D I G H I E RO , G . et al. (2020). Identiﬁcation of single spectral lines through supervised machine learning in a large HST survey (WISP): a pilot study for Euclid and WFIRST. The Astr ophysical Journal Supplement Series 249 12. https://doi.org/10.3847/1538- 4365/ab9a3a B A RO N C H E L L I , I . , S C A R L A TA , C . M . , R O D R I G U E Z - M U Ñ O Z , L . et al. (2021). Identiﬁcation of single spectral lines in large spectroscopic surveys using UMLA UT: an unsupervised machine learning algorithm based on unbiased topology. The Astr ophysical Journal Supplement Series 257 67. https://doi.org/10.3847/1538- 4365/ ac250c B AY E R , A . E . and S E L JA K , U . (2020). The look-elsewhere ef fect from a uniﬁed Bayesian and frequentist per- spectiv e. Journal of Cosmolo gy and Astr oparticle Physics 2020 009–009. https://doi.org/10.1088/1475- 7516/ 2020/10/009 B E R G E R , J . O . and P E R I C C H I , L . R . (1996). The intrinsic Bayes factor for model selection and prediction. Journal of the American Statistical Association 91 109–122. https://doi.org/10.2307/2291387 B O LTO N , A . S . , S C H L E G E L , D . J . , A U B O U R G , E . , B A I L E Y , S . et al. (2012). Spectral classiﬁcation and redshift measurement for the SDSS-III baryon oscillation spectroscopic survey. The Astr onomical Journal 144 144. https://doi.org/10.1088/0004- 6256/144/5/144 B Y R D , R . H . , L U , P . , N O C E DA L , J . and Z H U , C . (1995). A limited memory algorithm for bound constrained optimization. SIAM Journal on Scientiﬁc Computing 16 1190–1208. https://doi.org/10.1137/0916069 E U C L I D C O L L A B O R ATI O N (2025). Euclid. I. overvie w of the Euclid mission. Astr onomy & Astrophysics 697 A1. https://doi.org/10.1051/0004- 6361/202450810 D A V I E S , R . B . (1977). Hypothesis testing when a nuisance parameter is present only under the alternativ e. Biometrika 64 247. https://doi.org/10.2307/2335690 D A V I E S , R . B . (1987). Hypothesis testing when a nuisance parameter is present only under the alternativ e. Biometrika 74 33–43. https://doi.org/10.1093/biomet/74.1.33 G E Y E R , C . and T H O M P S O N , E . (1995). Annealing Markov chain Monte Carlo with applications to ancestral inference. Journal of the American Statistical Association 90 909–920. https://doi.org/10.1080/01621459. 1995.10476590 G R O S S , E . and V I T E L L S , O . (2010). Trial factors for the look else where effect in high energy physics. The Eur opean Physical Journal C 70 525–530. https://doi.org/10.1140/epjc/s10052- 010- 1470- 8 H A N S E N , B . E . (1996). Inference when a nuisance parameter is not identiﬁed under the null hypothesis. Econo- metrica 64 413–430. https://doi.org/10.2307/2171789 H A R R I S O N , I . , L O C H N E R , M . and B R OW N , M . L . (2017). Redshifts for galaxies in radio continuum surveys from Bayesian model ﬁtting of HI 21-cm lines. https://doi.org/10.48550/arXi v .1704.08278 H U B E RT Y , M . S . , N E D K OV A , K . V . , S ATTA R I , Z . et al. (2026). The NIRISS P ASSAGE spectroscopic redshift catalog in COSMOS. https://doi.org/10.48550/arXi v .2601.02220 22 J A M A L , S . , L E B RU N , V . , L E F È V R E , O . et al. (2018). Automated reliability assessment for spectroscopic red- shift measurements. Astr onomy and Astrophysics 611 A53. https://doi.org/10.1051/0004- 6361/201731305 L AU R E I J S , R . , A M I A U X , J . , A R D U I N I , S . et al. (2011). Euclid deﬁnition study report. ESA/SRE(2011)12 . arXiv:1110.3193. https://doi.org/10.48550/arXiv .1110.3193 M A R I N , J . - M . and R O B E R T , C . (2014). Bayesian Essentials with R , 2nd ed. Springer T exts in Statistics . Springer, New Y ork. P A R K , T., V A N D Y K , D . A . and S I E M I G I N O W S K A , A . (2008). Searching for narrow emission lines in X-ray spectra: computation and methods. The Astr ophysical Journal 688 807–825. https://doi.org/10.1086/591631 R E N C H E R , A . C . and S C H A A L J E , G . B . (2008). Linear Models in Statistics , 2nd ed. W iley-Interscience, Hobo- ken, N.J. R O B E RT S , G . O . and R O S E N T H A L , J . S . (1998). Optimal scaling of discrete approximations to Langevin diffusions. Journal of the Royal Statistical Society Series B: Statistical Methodology 60 255–268. https: //doi.org/10.1111/1467- 9868.00123 S Y E D , S . , B O U C H A R D - C Ô T É , A . , D E L I G I A N N I D I S , G . and D O U C E T , A . (2022). Non-reversible parallel tem- pering: a scalable highly parallel MCMC scheme. Journal of the Royal Statistical Society Series B: Statistical Methodology 84 321–350. https://doi.org/10.1111/rssb .12464 T O N RY , J . and D A V I S , M . (1979). A survey of galaxy redshifts. I - data reduction techniques. The Astronomical Journal 84 1511. https://doi.org/10.1086/112569 W I L K S , S . S . (1938). The large-sample distrib ution of the likelihood ratio for testing composite hypotheses. The Annals of Mathematical Statistics 9 60–62. https://doi.org/10.1214/aoms/1177732360 Z A M O R A , S . and D Í A Z , A . I . (2024). Re vising the cross correlation technique at high spectral resolution. https: //doi.org/10.48550/arXi v .2310.04133 B A YESIAN EMISSION LINE DETECTION 23 SUPPLEMENT TO “ A SCALABLE B A YESIAN FRAMEWORK FOR EMISSION LINE DETECTION AND REDSHIFT ESTIMA TION” S.1. Likelihood derivation. Recall that we impose orthogonality conditions on the lin- ear model in equation ( 7 ) in the main text. In particular , we require 1 ′ N [ S ( σ 2 )] − 1 X = 0 , 1 ′ N [ S ( σ 2 )] − 1 W sub ( ζ , δ ) = 0 , X ′ [ S ( σ 2 )] − 1 W sub ( ζ , δ ) = 0 , where the zero-v ector on the right-hand side of each equation above is of appropriate dimen- sion. T o do so, consider the projection matrix onto the span of 1 N as P 1 N = 1 N (1 ′ N [ S ( σ 2 )] − 1 1 N ) − 1 1 ′ N [ S ( σ 2 )] − 1 , and let M 1 N = I N − P 1 N denote the corresponding orthogonal complement projection ma- trix. Projecting X using this operator giv es X ⊥ = M 1 N X , which satisﬁes 1 ′ N [ S ( σ 2 )] − 1 X ⊥ = 0 . Next, deﬁne the projection matrix onto the span of X ⊥ as P X ⊥ = X ⊥ ( X ′ ⊥ [ S ( σ 2 )] − 1 X ⊥ ) − 1 X ′ ⊥ [ S ( σ 2 )] − 1 , and the corresponding orthogonal complement as M X ⊥ = I N − P X ⊥ . Then, for each ( ζ , δ ) , the orthogonalized matrix W sub , ⊥ ( ζ , δ ) = M X ⊥ M 1 N W sub ( ζ , δ ) satisﬁes the orthogonality conditions 1 ′ N [ S ( σ 2 )] − 1 W sub , ⊥ ( ζ , δ ) = 0 , X ′ ⊥ [ S ( σ 2 )] − 1 W sub , ⊥ ( ζ , δ ) = 0 . For simplicity of notation, we write X and W sub rather than X ⊥ and W sub , ⊥ , bearing in mind that this orthogonalization has been performed. W ith the orthogonality conditions above, the generalized least-squares estimators of α, β and η sub ( ζ ) for ﬁxed ( ζ , δ ) are ˆ α = (1 ′ N [ S ( σ 2 )] − 1 1 N ) − 1 1 ′ N [ S ( σ 2 )] − 1 y , (S.1) ˆ β = ( X ′ [ S ( σ 2 )] − 1 X ) − 1 X ′ [ S ( σ 2 )] − 1 ( y − ˆ α 1 N ) , (S.2) ˆ η sub ( ζ ) = [ W sub ( ζ , δ ) ′ [ S ( σ 2 )] − 1 W sub ( ζ , δ )] − 1 W sub ( ζ , δ ) ′ [ S ( σ 2 )] − 1 ( y − ˆ α 1 N − X ˆ β ) . (S.3) 24 Let r = y − ˆ α 1 N − X ˆ β . A routine calculation (see, e.g., Ch. 3 Marin and Robert , 2014 ) yields (S.4) p ( y | α, β , σ 2 , η sub ( ζ ) , ζ , δ ) ∝ exp ( − 1 2 ( η sub ( ζ ) − ˆ η sub ( ζ )) ′ × W sub ( ζ , δ ) ′ [ S ( σ 2 )] − 1 W sub ( ζ , δ )( η sub ( ζ ) − ˆ η sub ( ζ )) ) × exp  − 1 2 ( α − ˆ α )1 ′ N [ S ( σ 2 )] − 1 1 N ( α − ˆ α )  × exp  − 1 2 ( β − ˆ β ) ′ X ′ [ S ( σ 2 )] − 1 X ( β − ˆ β )  × exp ( − 1 2 [ r − W sub ( ζ , δ ) ˆ η sub ( ζ )] ′ × [ S ( σ 2 )] − 1 [ r − W sub ( ζ , δ ) ˆ η sub ( ζ )] ) × ( σ 2 ) − N/ 2 where the cross-terms are zero by the orthogonal construction and properties of least-squares estimates (and the prior on δ is a constant by uniformity). From the full model likelihood in equation ( S.4 ), we can obtain the lik elihood for the null model by setting η sub ( ζ ) = ˆ η sub ( ζ ) = 0 . S.2. Marginal density under alternativ e model. This section provides details for com- puting the density p ( y | ˆ σ 2 , ζ , δ ) described in Section 5.2.1 . T o do so, we consider the encom- passing model deﬁned in equation ( 13 ). S.2.1. Encompassing model derivation. In this setting, we ha ve the (trivial but ke y) re- lationship p ( η sub ( ζ ) | ζ ) = p ( η ∗ sub ( ζ ) | ζ )1 { η ∗ sub ( ζ ) ∈A ( ζ ) } Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ζ ) , with A ( ζ ) deﬁned as in ( 9 ) and with η ∗ sub ( ζ ) satisfying ( 13 ). Using the encompassing model, equation ( 15 ) giv es an expression for p ∗ ( y | ˆ σ 2 , ζ , δ ) after integrating out α, β , and η ∗ sub ( ζ ) directly due to conjugacy . Using this expression, we can integrate out α, β , and η sub ( ζ ) to obtain p ( y | ˆ σ 2 , ζ , δ ) . First, α and β can be inte grated out directly by conjugacy , so p ( y | ˆ σ 2 , η sub ( ζ ) , ζ , δ ) is just the density of a multi variate normal distrib ution. Next, η sub ( ζ ) can be integrated out as follo ws: p ( y | ˆ σ 2 , ζ , δ ) = Z p ( y | ˆ σ 2 , η sub ( ζ ) , ζ , δ ) p ( η sub ( ζ ) | ζ ) dη sub ( ζ ) = Z p ( y | ˆ σ 2 , η ∗ sub ( ζ ) , ζ , δ ) p ( η ∗ sub ( ζ ) | ζ )1 { η ∗ sub ( ζ ) ∈A ( ζ ) } Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ζ ) dη ∗ sub ( ζ ) = p ∗ ( y | ˆ σ 2 , ζ , δ ) Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ζ ) Z p ( η ∗ sub ( ζ ) | ˆ σ 2 , ζ , δ, y )1 { η ∗ sub ( ζ ) ∈A ( ζ ) } dη ∗ sub ( ζ ) = p ∗ ( y | ˆ σ 2 , ζ , δ ) Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ˆ σ 2 , ζ , δ, y ) Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ζ ) , B A YESIAN EMISSION LINE DETECTION 25 where p ∗ ( y | ˆ σ 2 , ζ , δ ) = Z p ( y | ˆ σ 2 , η ∗ sub ( ζ ) , ζ , δ ) p ( η ∗ sub ( ζ ) | ζ ) dη ∗ sub ( ζ ) . From the deriv ation in Appendix S.3.3 , it follo ws that η ∗ sub ( ζ ) | σ 2 , ζ , δ, y follows a multiv ari- ate normal distribution with known mean and cov ariance matrix. Next, we demonstrate how to quickly compute (S.5) Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | σ 2 , ζ , y ) and Pr( η ∗ sub ( ζ ) ∈ A ( ζ ) | ζ ) in R. S.2.2. Linear inequality constr aints. The constraints in the set A ( ζ ) can always be writ- ten as A ( ζ ) η ∗ sub ( ζ ) ≥ 0 for an appropriate matrix, A ( ζ ) . T o construct the latter, note that the nonnegati vity and ratio constraints in equation ( 9 ) are partly redundant. W e therefore retain only a minimal set of inequalities deﬁning the same feasible region for ﬁxed ζ . Whenever two or more signal intensities (the η k ’ s) are related by ratio constraints, keep only the non- negati vity constraint for the ﬁrst index in that ordered set (i.e. the smallest of k = 3 , 4 , 7 in K ( ζ ) ), and include only the pairwise ratio constraints comparing each subsequent element to the previous one. For example, when K ( ζ ) = { 3 , 4 , 5 , 6 , 7 } , the acti ve constraints can be written as η 4 ( ζ ) − 2 . 174 η 3 ( ζ ) ≥ 0 , η 7 ( ζ ) − 2 . 86 η 4 ( ζ ) ≥ 0 , η k ( ζ ) ≥ 0 , k ∈ K ( ζ ) . The reduced, nonredundant system is therefore η 4 ( ζ ) − 2 . 174 η 3 ( ζ ) ≥ 0 , η 7 ( ζ ) − 2 . 86 η 4 ( ζ ) ≥ 0 , η 3 ( ζ ) ≥ 0 , η 5 ( ζ ) ≥ 0 , η 6 ( ζ ) ≥ 0 , which we can express compactly in matrix form as A ( ζ ) η sub ( ζ ) ≥ 0 , where A ( ζ ) =         − 2 . 174 1 0 0 0 0 − 2 . 86 0 0 1 1 0 0 0 0 0 0 1 0 0 0 0 0 1 0         , η sub ( ζ ) =         η 3 η 4 η 5 η 6 η 7         . Since A ( ζ ) η ∗ sub ( ζ ) | ˆ σ 2 , ζ , y follo ws a multi variate normal distribution with mean A ( ζ ) m ˆ σ 2 , ˆ η sub ( ζ ) ,δ and cov ariance matrix A ( ζ ) V ˆ σ 2 , ˆ η sub ( ζ ) ,δ A ( ζ ) ′ , where m ˆ σ 2 , ˆ η sub ( ζ ) ,δ and V ˆ σ 2 , ˆ η sub ( ζ ) ,δ are the cen- trality and spread parameters for the truncated multiv ariate normal posterior of p ( η sub ( ζ ) | ˆ σ 2 , ζ , δ, y ) deﬁned in equations ( S.7 ) and ( S.8 ), respectiv ely . Importantly , these quantities are av ailable in closed form. The probability in ( S.5 ) corresponds to the nonnegati ve orthant probability for A ( ζ ) η ∗ sub ( ζ ) , which can be ev aluated directly using pmvnorm() in R. The same can be done using the prior mean and cov ariance instead. Like wise, from the prior in equation ( 13 ), we know that A ( ζ ) η ∗ sub ( ζ ) | ζ follo ws a multi variate normal distrib ution as well with mean A ( ζ ) a 2 ( ζ ) and cov ariance matrix A ( ζ ) B 2 ( ζ ) A ( ζ ) ′ , where a 2 ( ζ ) and B 2 ( ζ ) are as deﬁned in equation ( 14 ). S.3. Posterior expectations. In this section we deri ve the posterior mean for α, β , δ, and η sub ( ζ ) , conditional on σ 2 = ˆ σ 2 and ζ = ˆ ζ = arg max r p ( ζ r | y , ˆ σ 2 , M 1 ) . 26 S.3.1. Backgr ound components. From the likelihood deriv ation in Appendix S.2 , we write p ( α | ˆ σ 2 , ˆ ζ , y ) ∝ exp  − 1 2 ( α − ˆ α )1 ′ N [ S ( ˆ σ 2 )] − 1 1 N ( α − ˆ α )  p ( α ) ∝ exp  − 1 2 ( α − ˆ α )1 ′ N [ S ( ˆ σ 2 )] − 1 1 N ( α − ˆ α )  exp  − 1 2 b 2 0 ( α − a 0 ) 2  ∝ exp  − 1 2 s 2 α ( α − m α ) 2  , by completing the square, where ˆ α is deﬁned in equation ( S.1 ), and s 2 α =  1 ′ N [ S ( ˆ σ 2 )] − 1 1 N + 1 b 2 0  − 1 , m α = s 2 α  ˆ α 1 ′ N [ S ( ˆ σ 2 )] − 1 1 N + a 0 b 2 0  . In particular , α | ˆ σ 2 , ˆ ζ , y ≡ α | ˆ σ 2 , y ∼ N ( m α , s α ) . The deriv ation for β | σ 2 , ˆ ζ , y is similar . W e have p ( β | ˆ σ 2 , ˆ ζ , y ) ∝ exp  − 1 2 ( β − ˆ β ) ′ X ′ [ S ( ˆ σ 2 )] − 1 X ( β − ˆ β )  exp  − 1 2 b 2 1 β ′ β  ∝ exp n ( β − m β ) ′ S − 1 β ( β − m β ) o , where ˆ β is deﬁned in equation ( S.2 ), and S β =  X ′ [ S ( ˆ σ 2 )] − 1 X + 1 b 2 1 I J  − 1 , m β = S − 1 β  X ′ [ S ( ˆ σ 2 )] − 1 X ˆ β  , where I J is the J × J identity matrix. That is, β | ˆ σ 2 , ˆ ζ , δ, y ≡ β | ˆ σ 2 , y ∼ N J ( m β , S β ) . S.3.2. Emission line pr oﬁle width. The posterior expectation for δ is obtained as follows. For ﬁx ed ˆ σ 2 , ˆ ζ , the posterior probability of each δ ℓ is (S.6) p ( δ ℓ | ˆ σ 2 , ˆ ζ , y ) = p ( y | ˆ σ 2 , ˆ ζ , δ ℓ ) p ( δ ℓ )∆ δ P D ℓ =1 p ( y | ˆ σ 2 , ˆ ζ , δ ℓ ) p ( δ ℓ )∆ δ , where each p ( y | ˆ σ 2 , ˆ ζ , δ ℓ ) term is calculated according to equation ( 16 ). Therefore, the pos- terior conditional mean of δ is E ( δ | ˆ σ 2 , ˆ ζ , y ) = D X ℓ =1 δ ℓ p ( δ ℓ | ˆ σ 2 , ˆ ζ , y ) . S.3.3. Signal intensities. Next, we compute the conditional posterior mean, E ( η sub ( ˆ ζ ) | ˆ σ 2 , ˆ ζ , y ) . T o start, we deriv e the full conditional for η sub ( ζ ) when the prior on η sub | ζ is assumed to be the truncated multiv ariate normal deﬁned in equation ( 10 ). Using the likeli- hood in equation ( S.4 ), expressed in terms of the weighted least-squares estimate ˆ η sub ( ζ ) , the B A YESIAN EMISSION LINE DETECTION 27 conditional posterior density of η sub ( ζ ) giv en ˆ σ 2 , ζ , and δ can be written as p ( η sub ( ζ ) | ˆ σ 2 , ζ , δ, y ) ∝ exp ( − 1 2 ( η sub ( ζ ) − ˆ η sub ( ζ )) ′ × W sub ( ζ , δ ) ′ S ( ˆ σ 2 ) − 1 W sub ( ζ , δ )( η sub ( ζ ) − ˆ η sub ( ζ )) ) × exp  − 1 2 ( η sub ( ζ ) − a 2 ( ζ )) ′ B 2 ( ζ ) − 1 ( η sub ( ζ ) − a 2 ( ζ ))  × 1 { η sub ( ζ ) ∈A ( ζ ) } ∝ exp ( − 1 2 ( η sub ( ζ ) − m ˆ σ 2 , ˆ η sub ( ζ ) ,δ ) ′ V − 1 ˆ σ 2 , ˆ η sub ( ζ ) ,δ × ( η sub ( ζ ) − m ˆ σ 2 , ˆ η sub ( ζ ) ,δ ) ) 1 { η sub ( ζ ) ∈A ( ζ ) } , where m ˆ σ 2 , ˆ η sub ( ζ ) ,δ := V ˆ σ 2 , ˆ η sub ( ζ ) ,δ  W sub ( ζ , δ ) ′ [ S ( ˆ σ 2 )] − 1 W sub ( ζ , δ ) ˆ η sub ( ζ ) + B − 1 2 ( ζ ) a 2 ( ζ )  , (S.7) and (S.8) V ˆ σ 2 , ˆ η sub ( ζ ) ,δ :=  W sub ( ζ , δ ) ′ [ S ( ˆ σ 2 )] − 1 W sub ( ζ , δ ) + B − 1 2 ( ζ )  − 1 . Due to the restriction that η sub ( ζ ) ∈ A ( ζ ) , the posterior mean of η sub ( ˆ ζ ) | ˆ σ 2 , ˆ ζ , y is not simply m ˆ σ 2 , ˆ η sub ( ˆ ζ ) ,δ (like wise for V ˆ σ 2 , ˆ η sub ( ˆ ζ ) ,δ and the posterior variance). Instead, we consider the linear inequality constraints in equation ( 9 ), and need to compute (S.9) E [ η sub ( ˆ ζ ) | ˆ σ 2 , ˆ ζ , δ, y ] = E [ η ∗ sub ( ˆ ζ ) | A ( ˆ ζ ) η ∗ sub ( ˆ ζ ) ≥ 0 , ˆ σ 2 , ˆ ζ , δ, y ] , where η ∗ sub ( ˆ ζ ) denotes the corresponding non-truncated multi variate normal random variable in equation ( 13 ). T o simplify notation, we will suppress the conditioning on ˆ σ 2 , ˆ ζ , and y , since they are ﬁxed in what follows, and write T ( ˆ ζ ) = A ( ˆ ζ ) η ∗ sub ( ˆ ζ ) , ˆ m δ = m ˆ σ 2 , ˆ η sub ( ˆ ζ ) ,δ , and ˆ V δ = V ˆ σ 2 , ˆ η sub ( ˆ ζ ) ,δ . The right-hand side of ( S.9 ) can then be expressed as E [ η ∗ sub ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , δ ] . T o ev aluate this expectation, note that η ∗ sub ( ˆ ζ ) and T ( ˆ ζ ) are jointly multiv ariate normal (conditional on δ ) with mean vector ( ˆ m δ , A ( ˆ ζ ) ˆ m δ ) ′ and cov ariance matrix  ˆ V δ ˆ V δ A ( ˆ ζ ) ′ A ( ˆ ζ ) ˆ V δ A ( ˆ ζ ) ˆ V δ A ( ˆ ζ ) ′  . Let u ∈ R denote an arbitary ﬁxed v alue of T ( ˆ ζ ) . F or any such u , it can be sho wn (e.g., Theorem 4.4d, Rencher and Schaalje , 2008 ) that the conditional mean is E [ η ∗ sub ( ˆ ζ ) | T ( ˆ ζ ) = u, δ ] = ˆ m δ + ˆ V δ A ( ˆ ζ ) ′ [ A ( ˆ ζ ) ˆ V δ A ( ˆ ζ ) ′ ] − 1 [ u − A ( ˆ ζ ) ˆ m δ ] . By the law of iterated e xpectation, E [ η ∗ sub ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , δ ] = E { E [ η ∗ sub ( ˆ ζ ) | T ( ˆ ζ ) = u, δ ] | T ( ˆ ζ ) ≥ 0 } = ˆ V δ A ( ˆ ζ ) ′ [ A ( ˆ ζ ) ˆ V δ A ( ˆ ζ ) ′ ] − 1 [ E ( T ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , δ ) − A ( ˆ ζ ) ˆ m δ ] + ˆ m δ , 28 in which E ( T ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , δ ) is the mean of a multi variate normal distribution truncated to the positiv e orthant, and can be efﬁciently computed in R, using tmvtnorm::mtmvnorm . Next, since δ is treated as a discrete parameter, equation ( S.9 ) reduces to: E ( η sub ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , ˆ σ 2 , ˆ ζ , y ) = D X ℓ =1 E ( η ∗ sub ( ˆ ζ ) | T ( ˆ ζ ) ≥ 0 , ˆ σ 2 , ˆ ζ , δ ℓ , y ) p ( δ ℓ | ˆ σ 2 , ˆ ζ , y ) , using the posterior weights for δ ℓ deﬁned in equation ( S.6 ). S.4. V ariance estimation. W e obtain a plugin estimate of σ 2 by maximizing the likeli- hood in equation ( S.4 ). For the optimization done here, δ and ζ are not discretized since no integration is necessary . Additionally , only linear equality constraints are imposed on η sub ( ζ ) since the inequality constraints in equation ( 9 ) are imposed through the prior in our model. The estimation of σ 2 proceeds as follo ws. Initialize σ 2 , ζ , and δ at a point, ( σ 2 (0) , ζ (0) , δ (0) ) , where the subscript denotes the iteration index. For iteration t = 1 , 2 , . . . , repeat the following steps until con vergence: 1. Gi ven ( σ 2 ( t − 1) , ζ ( t − 1) , δ ( t − 1) ) , compute the generalized least-squares estimates of α, β , and η sub ( ζ ) deﬁned in equations ( S.1 )-( S.3 ). Denote these as ( ˆ α ( t ) , ˆ β ( t ) , ˆ η sub , ( t ) ( ζ ( t − 1) )) . 2. Holding ( ˆ α ( t ) , ˆ β ( t ) , ˆ η sub , ( t ) ( ζ ( t − 1) )) ﬁxed, maximize the likelihood in ( S.4 ) with respect to ( σ 2 , ζ , δ ) using numerical optimization yielding ( σ 2 ( t ) , ζ ( t ) , δ ( t ) ) . The optimization in Step 2 is performed using the L-BFGS-B method via optim in R, sub- ject to box constraints σ 2 ∈ (0 . 2 , ∞ ) , ζ ∈ (0 , 4] and δ ∈ [1 , 3] . Conv ergence is determined by the default stopping criteria of optim for L-BFGS-B ( Byrd et al. , 1995 ), which terminates when the relati ve change in the objectiv e function between successi ve iterations falls below approximately 10 − 9 , or stop when the iteration limit ( t = 100 ) is reached. The resulting esti- mate, ˆ σ 2 , is the ﬁnal iterate σ 2 ( t ) . Estimates of the remaining parameters obtained during this procedure are not used in subsequent analysis. S.5. Posterior predicti ve samples. W e describe here the procedure for generating pos- terior predicti ve samples for the model with likelihood and prior deﬁned in Sections 3 and 4 , respecti vely . Recall that σ 2 is assigned a point-mass empirical Bayes prior ﬁxed at σ 2 = ˆ σ 2 . The posterior can be factorized as p ( θ | y ) ∝ p ( η sub ( ζ ) | α , β , ˆ σ 2 , ζ , δ, y ) p ( α | ˆ σ 2 , y ) p ( β | ˆ σ 2 , y ) p ( ζ , δ | ˆ σ 2 , y ) , where the conditional distributions of α and β are normal, and η sub | α, β , ˆ σ 2 , ζ , δ, y follows a truncated multi variate normal distrib ution (see Appendix S.3 ). W e sample ( ζ , δ ) o ver their discrete joint grid using the approximation ˆ p ( ζ r , δ ℓ | ˆ σ 2 , y ) = p ( y | ˆ σ 2 , ζ r , δ ℓ ) p ( δ ℓ ) p ( ζ r )∆ δ ∆ ζ P R r =1 P D ℓ =1 p ( y | ˆ σ 2 , ζ r , δ ℓ ) p ( δ ℓ ) p ( ζ r )∆ δ ∆ ζ , which is computed in the process of obtaining p ( ζ | ˆ σ 2 , y ) . A single posterior draw , θ ( t ) , is obtained sequentially: 1. Sample ( ζ ( t ) , δ ( t ) ) ∼ ˆ p ( ζ , δ | ˆ σ 2 , y ) , 2. Sample α ( t ) ∼ p ( α | ˆ σ 2 , y ) and β ( t ) ∼ p ( β | ˆ σ 2 , y ) (note these do not depend on the sam- ples obtained in Step 1). 3. Sample η ( t ) sub ( ζ ( t ) ) ∼ p ( η | α ( t ) , β ( t ) , ˆ σ 2 , ζ ( t ) , δ ( t ) , y ) using the function rtmvnorm in the tmvtnorm package in R with moments deﬁned above, and constraints imposed through A ( ζ ( t ) ) . The remaining components of η ( t ) ( ζ ( t ) ) are then obtained deterministically from the model constraints. B A YESIAN EMISSION LINE DETECTION 29 Finally , a posterior predicti ve draw , ˜ y ( t ) , is obtained by sampling from the likelihood con- ditioned on the posterior draw . That is, ˜ Y ( t ) ∼ p ( y | θ ( t ) ) . Repeating this process for t = 1 , ..., M yields a collection of posterior predictiv e samples suitable for model checking. S.5.1. P osterior r esiduals. F or model assessment, we examine the (standardized) resid- ual curves, deﬁned as wa velength-speciﬁc (standardized) residuals obtained by subtracting the posterior predictiv e mean from the observed data. For a single spectrum, at wav elength index i , deﬁne the standardized residual as (S.10) r i = y i − 1 M P M t =1 E ( ˜ Y ( t ) i | θ ( t ) , y ) p S ii ( ˆ σ 2 ) , where θ ( t ) denotes the t -th posterior sample, ˜ Y ( t ) i is the i -th entry of ˜ Y ( t ) , the posterior pre- dicti ve random vector corresponding to the t -th sample, E ( ˜ Y ( t ) i | θ ( t ) , y ) denotes the posterior predicti ve mean ﬂux value at the i -th wa velength, S ii ( ˆ σ 2 ) denotes the i -th diagonal entry of S ( ˆ σ 2 ) , and M is the total the number of posterior samples. The numerator therefore subtracts the posterior predictiv e mean (approximated by Monte Carlo averaging) from the observed ﬂux. One could alternatively use posterior predicti ve draws, ˜ y ( t ) j , t = 1 , . . . , M , directly in place of the conditional expectations. Howe ver , doing so would introduce additional Monte Carlo noise into the residuals and obscure assessment of the posterior predicti ve mean ﬁt. S.5.2. P osterior pr edictive p-values. For model assessment, we also consider posterior predicti ve p-values (PPPs). The speciﬁc PPPs used in Appendix S.6 are described here. W e distinguish between global PPPs, which produce a single scalar summary for each spectrum, and local PPPs which produce a vector of v alues indexed by w av elength. W e ﬁrst recall the deﬁnition of a global PPP . Gi ven posterior draws θ ( t ) , for t = 1 , . . . , M , and the associated posterior predicti ve samples, ˜ y ( t ) , let T ( y , θ ) denote a discrepancy func- tion which is some measure of model ﬁt (e.g., mean squared error), which may depend on both the data and the parameters. The global PPP is then estimated via Monte Carlo as [ P P P = 1 T T X t =1 I { T ( ˜ y ( t ) ,θ ( t ) ) ≥ T ( y,θ ( t ) ) } , where I {·} denotes the indicator function, and y denotes the observ ed ﬂux vector . W e consider three global discrepancies: mean squared error ( T MSE ), the maximum standardized residual ( T Max ), and the standardized median absolute de viation from the median ( T MAD ). Formally , they are deﬁned as T MSE ( y , θ ) = 1 N N X i =1 ( y i − E ( ˜ Y i | θ , y )) 2 S ii ( ˆ σ 2 ) , (S.11) T Max ( y , θ ) = max i ( | y i − E ( ˜ Y i | θ , y ) | p S ii ( ˆ σ 2 ) ) , (S.12) T MAD ( y , θ ) = med i ( | y i − med ( ˜ Y i | θ , y ) | p S j j ( ˆ σ 2 ) ) , (S.13) where E ( ˜ Y i | θ , y ) and med ( ˜ Y i | θ , y ) represent the posterior predictiv e mean and median ﬂux v alue at the i -th wav elength, respectively . Next, we consider a local PPP based on squared 30 null alternative 1.00 1.25 1.50 1.75 2.00 2.25 −2 0 2 −2 0 2 Observed wav elength ( µ m) Standardized residual F I G S . 1 . Aggre gated standardized r esidual curves for the 213 labeled spectra, separated accor ding to log poste- rior odds model selection into 203 alternative cases (top) and 10 null cases (bottom). The dark blue curve shows the pointwise median of the standar dized r esiduals acr oss spectra, and the light blue band r epr esents the 2.5% and 97.5% pointwise quantiles. The greater variability in the null panel r eﬂects the smaller number of spectra classiﬁed as null. standardized residuals, where the j -th local PPP is gi ven by (S.14) T local ,i ( y , θ ) = ( y i − E ( ˜ Y i | θ , y )) 2 S ii ( ˆ σ 2 ) . That is, we simply compute the i -th entry in the summation deﬁning T MSE ( y , θ ) , but do not av erage across wav elength. These local PPPs are useful in that they allow us to identify indices i for which the local PPPs are small. For example, if i = 4 , 5 , 6 , 7 , 8 all lead to small PPPs, then we expect some structural model misﬁt for the corresponding w av elengths. S.6. Model checking. This appendix presents graphical and numeric diagnostics to as- sess model ﬁt. These includes residual curv es based on posterior samples posterior predicti ve p-v alue summaries (as deﬁned in Appendix S.5 ). Model checking for spectral data is some- what challenging duet to the presence of sharp emission features, which generate extreme observ ations. Residual-based diagnostics must therefore be interpreted with care, as large residuals may reﬂect emission lines that are slightly underﬁt rather than signiﬁcant model misﬁt overall. Below , we describe diagnostics designed to characterize model ﬁt in this set- ting. S.6.1. Residual curves: 213 labeled spectra. W e ﬁrst assess the mean ﬁt (background and signal components) across the 213 labeled spectra. For each spectrum, the per- wa velength residuals r i , i = 1 , . . . , N , deﬁned in Appendix S.5.1 , are index ed by wav elength and therefore form a standardized residual curve ov er the N observed wav elength values. Model ﬁt is ev aluated after model selection. Speciﬁcally , for each spectrum we select either the null model, M 0 , or the alternativ e model, M 1 , based on the log posterior odds. Residual B A YESIAN EMISSION LINE DETECTION 31 null alternative 1.00 1.25 1.50 1.75 2.00 2.25 −2 0 2 4 −2 0 2 4 Observed wav elength ( µ m) Standardized residual F I G S . 2 . Aggr egated standardized r esidual curves for the 2800 unlabeled spectra, separated accor ding to log posterior odds model selection into 937 alternative cases (top) and 1863 null cases (bottom). The dark blue curve shows the pointwise median of the standardized residuals across spectra, and the light blue band repr esents the 2.5% and 97.5% pointwise quantiles. The gr eater variability in the null panel reﬂects the smaller number of spectra classiﬁed as null. curves are then aggre gated separately within each selected model. Figure S.1 displays the ag- gregated curves as functions of wa velength across the 203 spectra classiﬁed as M 1 and the 10 spectra classiﬁed as M 0 . The posterior predicti ve means in equation ( S.10 ) are approximated using M = 1000 posterior draws for each spectrum. For both the null and alternativ e cases, ideal model ﬁt would produce residual curves whose median ﬂuctuates around zero, with the 95% pointwise quantile bands largely con- tained within ± 2 . Because the aggreg ation is based on only 203 spectra for M 1 and 10 spec- tra for M 0 , some variability around zero is expected. The null panel exhibits substantially more v ariability due to the small sample size. The primary diagnostic concern is systematic deviation from zero. Such behavior is most e vident near the boundaries of the observed wav elength range. In particular, near the edges and detector gaps the median residual curve shows noticeable de viations from zero, espe- cially under the alternati ve model. This reﬂects the limited information av ailable to constrain the background and emission line model components in these regions due to being near a boundary . A way from detector gaps, the residual curves generally beha ve like constant func- tions centered around zero with some random ﬂuctuation. S.6.2. Residual curves: 2800 unlabeled spectra. Figure S.2 sho ws aggregated residual curves for the 2800 unlabeled spectra, separated according to log posterior odds model selec- tion into 937 alternati ve cases and 1863 null cases. Because this dataset contains substantially more spectra classiﬁed under both models, the aggregated residual curves are smoother and more stable than those observed for the labeled case. For spectra with detected emission lines, similar edge ef fects are observed near detector gaps, consistent with the beha vior seen in the labeled data. These boundary ef fects are considerably less pronounced in the null case. 32 Null model Alternative model MAD Max MSE 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0 20 40 60 0 10 20 30 40 50 0 10 20 30 P oster ior predictive p−v alue Count F I G S . 3 . Histograms of posterior predictive p-values for the 213 labeled spectra, computed using the T MAD , T Max , and T MSE discr epancies in ( S.11 ) – ( S.13 ) (top to bottom). Spectra are ﬁrst classiﬁed accor ding to the log posterior odds model selection criterion, and goodness-of-ﬁt is then evaluated within the selected model. The left column corresponds to spectra classiﬁed as null and the right column to those classiﬁed as alternative. The null histograms ar e sparse due to the small number of spectra selected as null in the labeled set. In both panels, the median residual curves remain close to zero across most wav elengths, indicating that, after model selection, the selected model captures the overall mean structure of the spectra reasonably well within each group. Patterns in residual variability are e vident, ho wev er . Most notably , over the wa velength range [1 . 03 µ m , 1 . 27 µ m] , the spread of the stan- dardized residuals decreases monotonically in both the null and alternati ve cases, a pattern not observed else where. This suggests residual heteroskedasticity in the observed data be- yond what is captured by the reported measurement errors. Because the v ariance model takes the form σ 2 I N + C as in equation ( 6 ), heteroskedasticity enters only through the reported measurement errors, while the estimated component σ 2 is global. Consequently , wa velength-dependent miscalibration of the measurement errors is not corrected by the model and persists in the posterior predictiv e distribution. The contraction of the standardized residual bands ov er [1 . 03 µ m , 1 . 27 µ m] indicates that v ariability may be slightly overestimated there. This form of misspeciﬁcation is less concerning than systematic underestimation of v ariance, as it generally produces more conservati ve inference, including weaker e vidence for line detection and wider redshift HPD sets. S.6.3. P osterior predictive c hecks: 213 labeled spectra. Figure S.3 sho ws histograms of global posterior predictiv e p-values (PPPs) for the 213 labeled spectra. A substantial portion of spectra e xhibit extreme PPP values (near 0 or 1), which at ﬁrst glance suggests poor model ov erall ﬁt. Howe ver , these discrepancies are based on global summaries and are sensitiv e to isolated outliers and the quality of the v ariance estimation. In particular , PPPs near zero typically arise when the posterior predictiv e distribution f ails to cov er the observed ﬂux at one or more wa velengths. This commonly occurs in cases of B A YESIAN EMISSION LINE DETECTION 33 −5 0 5 10 1.00 1.25 1.50 1.75 2.00 2.25 Flux 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 Obser ved w avelength ( µ m) PPP F I G S . 4 . The top panel shows the observed data for spectrum 00004 (black dots) together with the pointwise median and 2.5% and 97.5% posterior predictive quantiles (dark and light blue). The bottom panel shows the corr esponding pointwise PPPs based on T local in Appendix S.5.1 , with dashed red lines at 0.05 and 0.95. A few wavelengths have PPPs below 0.05, but none ar e consecutive. 25 50 75 100 125 1.00 1.25 1.50 1.75 2.00 2.25 Flux 0.00 0.25 0.50 0.75 1.00 1.25 1.50 1.75 2.00 2.25 Obser ved w avelength ( µ m) PPP F I G S . 5 . The top panel shows the observed data for spectrum 00201 (black dots) together with the pointwise median and 2.5% and 97.5% posterior pr edictive quantiles (dark and light blue). The pr ominent bump is not captur ed by the posterior predictive intervals. The bottom panel shows the corresponding pointwise PPPs based on T local in ( S.14 ) , with dashed r ed lines at 0.05 and 0.95. A sequence of PPPs below 0.05 aligns with this lack of ﬁt, and the PPPs ar e shifted above 0.5, indicating inﬂated posterior pr edictive variance at each wavelength. slight underﬁtting of strong emission lines. Because the global discrepancies aggre gate across wa velengths, a single large residual can driv e the PPP close to zero. For example, spectrum 34 Null model Alternative model MAD Max MSE 0.00 0.25 0.50 0.75 1.00 0.00 0.25 0.50 0.75 1.00 0 50 100 150 0 50 100 150 200 250 0 100 200 P oster ior predictive p−v alue Count F I G S . 6 . Histograms of posterior predictive p-values for the 2800 unlabeled spectra, computed using the T MAD , T Max , and T MSE discr epancies in ( S.11 ) – ( S.13 ) (top to bottom). Spectra are ﬁrst classiﬁed accor ding to the log posterior odds model selection criterion, and goodness-of-ﬁt is then evaluated within the selected model. The left column corr esponds to spectra classiﬁed as null and the right column to those classiﬁed as alternative . 00033 in Figure 6 of the main text has a T MSE -based PPP of zero, ev en though the o ver - all posterior predictiv e ﬁt is satisfactory . Conv ersely , PPP values near one indicate that the predicti ve discrepancies are systematically larger than the observ ed discrepancy . This behav- ior is often consistent with predictiv e o verdispersion due, for e xample, to inﬂated variance estimates. In this setting, such misspeciﬁcation is less concerning, as it primarily leads to conserv ativ e inference as noted previously . The presence of strong emission lines itself can inﬂate v ariance estimates, making this behavior unsurprising. Because the global PPPs are sensiti ve to isolated de viations, we focus on structural model failures using local PPPs based on ( S.14 ) computed at each wa velength. These allo w identiﬁ- cation of systematic discrepancies across contiguous wa velength regions, such as underﬁtting of emission lines or background misspeciﬁcation. T o ﬂag such behavior , we identify spec- tra exhibiting multiple consecutive local PPPs near zero (e.g., we ﬂag any spectra with four consecuti ve PPPs belo w 0.05). The top panel of Figure S.4 shows the posterior predictiv e distribution for spectrum 00004, which appears to ﬁt well according to visual assessment as discussed in the main text. The bottom panel of the ﬁgure sho ws the corresponding sequence of local PPPs across wa ve- length, with dashed red lines marking extreme range cutoffs of 0.05 and 0.95 (chosen some- what arbitrarily). Although a few wav elengths yield small PPP values, they are isolated rather than consecutiv e. The corresponding global PPPs ( T MAD , T Max , and T MSE ) are 0.97, 0.59, and 1, respecti vely , for reference. In contrast, Figure S.5 shows the same two plots for spectrum 00201, which exhibits a prominent emission line feature not captured by the posterior predictiv e distribution. This corresponds to a sequence of consecutive small local PPPs aligned with the structural dis- crepancy . The global PPPs for this spectrum are 1, 0, and 0.84 for T MAD , T Max , and T MSE , respecti vely , reﬂecting both localized underﬁtting and v ariance inﬂation. B A YESIAN EMISSION LINE DETECTION 35 1.0 1.2 1.4 1.6 1.8 2.0 2.2 −10 10 30 50 Spectrum 00057 1.0 1.2 1.4 1.6 1.8 2.0 2.2 0 500 1000 Spectrum 00230 1.0 1.2 1.4 1.6 1.8 2.0 2.2 4000 6000 8000 Spectrum 02802 1.0 1.2 1.4 1.6 1.8 2.0 2.2 1500 2500 3500 Spectrum 03006 Obser ved w avelength ( µ m) Flux F I G S . 7 . Plots of the observed ﬂux values versus wavelength for four low quality spectra out of the unlabeled set of 2800. All four of these spectra ar e ﬂag ged as having poor model ﬁt, and ar e r emoved fr om the main r esults due to obvious err ors in the detection and data recor ding pr ocess. T aken together , these examples illustrate that extreme global PPP values often reﬂect local- ized deviations rather than widespread model f ailure. Local PPPs provide a more informati ve diagnostic by identifying contiguous regions of misﬁt. In total, for the 213 labeled spectra, using a identiﬁcation rule of ﬁv e consecutiv e local PPPs less than or equal to 0.05, there are nine spectra identiﬁed for follow-up due to potentially poor model ﬁt. Note that while the posterior predicti ve ﬁt is poor , across the nine ﬂagged spectra, the largest absolute difference between the marginal MAP redshift estimate and the reported redshift from astronomers is 0.00187. That is, these ﬁts do not lead to pathological redshift ﬁtting as they still identify the locations of emission lines e ven if heights of emission lines are underﬁt, or the variance is poorly estimated. S.6.4. P osterior pr edictive checks: 2800 unlabeled spectra. Figure S.6 sho ws histograms of global posterior predicti ve p-v alues (PPPs) for the 2800 unlabeled spectra. For spectra se- lected as containing emission lines, the PPP behavior closely mirrors that observed in the labeled sample. Under the null model, ho wev er , the substantially lar ger number of spectra al- lo ws a clearer assessment of discrepancy beha vior . The T MSE and T MSE discrepancies remain sensiti ve to isolated extreme residuals (as they should), often producing PPPs near 0 or 1. In contrast, the T MAD discrepancy is more robust to outliers. Under the null model, the T MAD - based PPPs are concentrated around 0.5 with an approximately symmetric decline to ward 0 and 1, consistent with PPP behavior under “adequate” model ﬁt. Posterior predicti ve p-v alues are kno wn to exhibit conserv atism and do not follow a uni- form distribution under the null. Nevertheless, the concentration of T MAD -based PPPs near 0.5 in the null case suggests that the background and noise components of the model are gen- erally well calibrated when no emission lines are present. The more extreme values observed 36 for the T MSE and T Max discrepancies primarily reﬂect sensitivity to localized extremes in the observed spectra. Using a local PPP summary instead, we ﬂag spectra that contain at least one sequence of ﬁ ve consecutiv e local PPPs below 0.05. This threshold is some what arbitrary , b ut we saw in the high-quality data setting that it meaningfully ﬂagged poor structural model ﬁt. This criterion identiﬁes 85 of the 2800 spectra as exhibiting poor model ﬁt. Figure S.7 displays four representati ve e xamples. The ﬂagged spectra exhibit se veral recurring data quality issues. Spectrum 00057 contains an e xtended stretch of ﬂux values ﬁxed at zero, likely reﬂecting truncation or missing val- ues. Spectrum 00230 shows an abrupt jump in ﬂux v alues, increasing from values between 0 and 100 to values near 1500, consistent with detector or data-processing artifacts. Spectra 02802 and 03006 contain extremely lar ge ﬂux values (exceeding 8000 and 3000, respec- ti vely), likely caused by contamination from a nearby bright source rather than the tar get galaxy . These examples illustrate structural data errors that would likely compromise inference un- der any modeling procedure. W e therefore remove all 85 ﬂagged spectra from the unlabeled analysis. While this ﬁltering rule may exclude a small number of spectra with reliable data but imperfect model ﬁt, in general the observed misﬁt appears to be dri ven primarily by data corruption. In this setting, removing suspect spectra is preferable to retaining observations that would yield unreliable redshift estimates. S.7. Prior h yperparameter sensitivity analysis. In the sensiti vity analysis that follows, each subsection examines variation in the hyperparameters of a single model component, while holding all other prior hyperparameters ﬁxed at the values described in Section 4 and used in the main analysis. For example, when assessing sensitivity to the signal intensity prior , the intercept prior , background speciﬁcation (including the number of components and prior variances), line width prior, and redshift prior are held ﬁxed at their values used in the main analysis. Sensiti vity results are conducted using only the 213 high-quality spectra. For each prior hyperparameter setting, we report the following summaries: (i) histograms of the log posterior odds; (ii) a table of rejection rates based on the log posterior odds, and false positiv e rates on simulated data; (iii) histograms of the resulting redshift mar ginal MAP estimates; and (i v) boxplots of total HPD widths. T o estimate false positiv e rates (FPRs) under each hyperparameter conﬁguration, we gen- erate synthetic null datasets as follows. For each of the 213 labeled spectra, we ﬁrst ﬁt the proposed model under the null hypothesis using the baseline priors from the main analy- sis. Using the resulting posterior estimates under the null model, we generate one synthetic dataset from the corresponding null-model likelihood. This produces 213 synthetic null spec- tra, one for each original JWST spectrum. The synthetic null data are generated using the baseline prior speciﬁcation rather than re- ﬁtting under each alternativ e hyperparameter conﬁguration. This ensures that the null data- generating mechanism remains ﬁxed across sensiti vity settings. For each prior hyperparam- eter conﬁguration under consideration, we then reapply the full detection procedure to the 213 synthetic null spectra and record the proportion of cases in which the alternative model is selected. This proportion serves as an empirical estimate of the FPR under that prior spec- iﬁcation. S.7.1. Alpha hyperparameter . The results in the main paper were replicated with the priors α ∼ N (5 , 10 2 ) and α ∼ N (0 , 100 2 ) . The results were numerically identical to the main analysis and are thus omitted. B A YESIAN EMISSION LINE DETECTION 37 Setting 4 Setting 5 Setting 6 Setting 1 Setting 2 Setting 3 0 100 200 300 0 100 200 300 0 100 200 300 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 0 10 20 30 40 0 20 40 0 20 40 Log posterior odds Count F I G S . 8 . Each histogr am depicts the log posterior odds for a differ ent hyperparameter setting for the number of backgr ound basis functions and beta coefﬁcient prior variances. A small number of larg e log posterior odds values ar e cutoff fr om each plot. The dashed blue line at zer o indicates the cutoff thr eshold. Setting Rejection rate (Signal) FPR (Null) 1 0.873 0.019 2 0.962 0.019 3 0.981 0.901 4 0.873 0.019 5 0.962 0.019 6 0.981 0.901 T A B L E S . 1 Rejection pr obabilities on JWST data labeled as containing a signal by astronomer s and the false positive rate (FPR) on simulated null data. Sample size of n = 213 per setting; maximum Monte Carlo standard err or ≤ 0 . 023 . S.7.2. Beta hyperparameter and number of bases. T o assess the sensitivity to the ﬂex- ibility of the background function, we v ary both the number of basis functions, J , and the prior v ariance for each β j . Speciﬁcally , we consider the following cases: 1. J = 30 with b 2 1 = 30 2 2. J = 15 with b 2 1 = 30 2 3. J = 5 with b 2 1 = 30 2 4. J = 30 with b 2 1 = 10 2 5. J = 15 with b 2 1 = 10 2 6. J = 5 with b 2 1 = 10 2 . Setting 2 ( J = 15 , b 1 = 30 ) corresponds to the prior speciﬁcation used in the main analysis. S.7.2.1. Impact on log posterior odds and FPRs. Figure S.8 displays six histograms of the log posterior odds corresponding to the hyperparameter settings described abov e. In the 38 Setting 4 Setting 5 Setting 6 Setting 1 Setting 2 Setting 3 0 1 2 3 4 0 1 2 3 4 0 1 2 3 4 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 Redshift Count F I G S . 9 . Each histogr am depicts the distribution of mar ginal MAP r edshift estimates for a differ ent hyperparam- eter setting for the number of backgr ound basis functions and beta coefﬁcient prior variances. top ro w of the ﬁgure, mo ving from left to right (with ﬁx ed prior variance b 2 1 = 30 2 ), decreas- ing the number of basis functions, J , shifts the distrib ution of the log posterior odds to wards positi ve values. This shift is practically meaningful, as it results in a larger number of spectra being identiﬁed as containing emission lines. The bottom row e xhibits an identical pattern for b 2 1 = 10 2 . That is, changing the prior variance from 30 2 to 10 2 has no ef fect on the log posterior odds sho wn here. T able S.1 reports the rejection rates for the 213 JWST spectra and the corresponding false positi ve rates (FPRs) based on 213 simulated null spectra. The results again sho w that sensi- ti vity to the background speciﬁcation is dictated by the number of basis functions J . When J = 5 basis functions are used, the FPR increases dramatically to 0.901, indicating se vere anti-conserv ativ e behavior arising from an inﬂexible background model. The accompanying increase in rejection rate is therefore misleading. In contrast, for J = 15 and J = 30 , the FPR remains low and stable, and is unaf fected by the choice of prior variance in this range. Al- though the FPR is identical for J = 15 and J = 30 , the J = 30 conﬁguration yields slightly lo wer rejection rates ov erall, showing that the detection power is mildly sensiti ve to the de- gree of background ﬂexibility (as w as observed with the log posterior odds distrib utions). S.7.2.2. Impact on r edshift distributions. Next, we examine the sensiti vity of redshift estimation to changes in background hyperparameters. Figure S.9 presents six histograms of marginal MAP redshift estimates for all 213 spectra. In contrast to the main analysis, we do not restrict attention to spectra identiﬁed as containing emission lines, in order to assess the ov erall impact of the hyperparameter settings on redshift v ariability . Once again, the results are effecti vely identical across the prior variance choices, and depend primarily on J . Across J = 5 , 10 , 15 , the redshift distributions are broadly similar , though some dif ferences are visible. For example, the J = 30 conﬁguration exhibits a slightly higher concentration of lo w-redshift estimates relativ e to the other cases. B A YESIAN EMISSION LINE DETECTION 39 0.001 0.010 0.100 1.000 Setting 1 Setting 2 Setting 3 Setting 4 Setting 5 Setting 6 Hyperparameter combination T otal HPD width F I G S . 1 0 . Each histogr am depicts the total HPD width (on the log scale) for a differ ent hyperparameter setting for the number of backgr ound basis functions and beta coefﬁcient prior variances. Case IIIa Case IV a Case Ia Case IIa 0 100 200 300 0 100 200 300 0 5 10 15 20 0 10 20 30 0 5 10 15 20 0 5 10 15 20 Log posterior odds Count F I G S . 1 1 . Each histogram depicts the log posterior odds for a differ ent hyperparameter setting for the signal intensities. A small number of larg e log posterior odds values ar e cutoff from each plot. The dashed blue line at zer o indicates the cutoff thr eshold. S.7.2.3. Impact on r edshift HPD set widths. Finally , Figure S.10 shows boxplots of the total 99.865% HPD widths for redshift under each hyperparameter conﬁguration (displayed on the log scale). As before, results are identical across prior v ariance settings and depend only on J . The median HPD width is nearly identical for J = 15 and J = 30 , with a modest 40 Setting Rejection rate (Signal) FPR (Null) Ia 0.962 0.019 IIa 0.915 0.005 IIIa 0.915 0.005 IV a 0.977 0.164 T A B L E S . 2 Rejection pr obabilities on JWST data labeled as containing a signal by astronomer s and the false positive rate (FPR) on simulated null data. Sample size of n = 213 per setting; maximum Monte Carlo standard err or ≤ 0 . 025 . increase when J = 5 . Overall, while the HPD widths exhibit some variation across hyperpa- rameter settings, they are less sensiti ve than the false positi ve rates. S.7.3. Signal intensity hyperparameters. W e next assess the sensiti vity of the results to the signal intensity h yperparameters, which prove to be the most inﬂuential component of the prior speciﬁcation. W e consider the following cases: • Case Ia - Science informed prior (used in main text): a 2 = (0 . 02 , 0 . 02 , 0 . 0092 , 0 . 02 , 0 . 2 , 0 . 596 , 0 . 2 , 0 . 02 , 0 . 02 , 0 . 0271 , 0 . 02 , 0 . 05 , 0 . 02) ′ , b 2 = (0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 1 , 1 , 1 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01) ′ • Case IIa - Prior with scientiﬁc informed mean and large equal standard de viations: a 2 = (0 . 02 , 0 . 02 , 0 . 0092 , 0 . 02 , 0 . 2 , 0 . 596 , 0 . 2 , 0 . 02 , 0 . 02 , 0 . 0271 , 0 . 02 , 0 . 05 , 0 . 02) ′ , b 2 = (1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1) ′ • Case IIIa - Priors with equal means and large equal standard de viations: a 2 = (0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02) ′ , b 2 = (1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1 , 1) ′ • Case IV a - Priors with equal means and small equal standard de viations: a 2 = (0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02 , 0 . 02) ′ , b 2 = (0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01 , 0 . 01) ′ S.7.3.1. Impact on log posterior odds and FPRs. Figure S.11 shows histograms of the log posterior odds under the four prior settings described in the previous section. Cases IIa and IIIa are noticeably more conservati ve than the science-informed prior (Case Ia). Both impose larger and equal prior variances across emission lines, increasing ﬂexibility in the alternati ve model and thereby introducing a lar ger penalty in the mar ginal density of y under the alternativ e model. As a result, the log posterior odds shift downward relativ e to Case Ia. Case IV a, in contrast, is the most liberal conﬁguration, identifying only ﬁve spectra as null cases. Ho wev er, T able S.2 shows that the increased detection rate comes at a cost: the FPR is 0.164, compared to 0.019, 0.005, and 0.005 for Cases Ia-IIIa, respectiv ely . The tight prior constraints in Case IV a reduce the effecti ve penalty for the alternative model, making spu- rious detections more likely . Overall, while many spectra yield consistent decisions across prior choices, detection rates vary meaningfully across conﬁgurations. Increasing prior vari- ance generally leads to more conserv ati ve behavior , whereas ov erly tight priors inﬂate false positi ves. It is reassuring to see that the science-informed prior , chosen a priori based on expert belief, balances this trade-of f. B A YESIAN EMISSION LINE DETECTION 41 Case IIIa Case IV a Case Ia Case IIa 0 1 2 3 4 0 1 2 3 4 0 5 10 15 0 10 20 0 5 10 15 20 0 5 10 15 Redshift Count F I G S . 1 2 . Each histogram depicts the distribution of mar ginal MAP redshift estimates for a differ ent hyperpa- rameter setting for the signal intensities. 0.001 0.010 0.100 1.000 Case Ia Case IIa Case IIIa Case IV a Hyper parameter combination T otal HPD width F I G S . 1 3 . Each histogr am depicts the total HPD width (on the log scale) for a differ ent hyperparameter setting for the signal intensities. S.7.3.2. Impact on r edshift distributions. Figure S.12 shows histograms of the marginal MAP redshift estimates under each prior speciﬁcation, for all 213 spectra, including non- detections. The signal intensity prior has a pronounced ef fect on the o verall redshift distri- bution. Cases Ia and IV a produce broadly similar patterns, though Case IV a yields fewer redshifts near zero and four . In contrast, Cases IIa and IIIa exhibit a substantial decrease in 42 Case IIIb Case IVb Case Ib Case IIb 0 100 200 300 0 100 200 300 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 0 5 10 15 20 Log posterior odds Count F I G S . 1 4 . Each histogram depicts the log posterior odds for a differ ent hyperparameter setting for redshift pa- rameter . A small number of lar ge log posterior odds values ar e cutoff fr om each plot. The dashed blue line at zer o indicates the cutoff thr eshold. Setting Rejection rate (Signal) FPR (Null) Ib 0.962 0.019 IIb 0.977 0.033 IIIb 0.953 0.009 IVb 0.962 0.047 T A B L E S . 3 Rejection pr obabilities on JWST data labeled as containing a signal by astronomer s and the false positive rate (FPR) on simulated null data. Sample size of n = 213 per setting; maximum Monte Carlo standard err or ≤ 0 . 015 . redshifts near two, with man y more redshifts at lower and higher redshifts. This sensiti vity is practically important. Relaxing prior information about relativ e emission line strengths alters which line conﬁgurations are fa vored, directly af fecting the inferred redshift modes. S.7.3.3. Impact on r edshift HPD set widths. Figure S.13 shows boxplots for the total 99.865% HPD set widths (on the log scale). Cases Ia and IV a again exhibit similar behavior , with Case IV a showing a modest increase in median width. Cases IIa and IIIa display larger median HPD set widths overall, indicating greater posterior uncertainty in redshift. Interest- ingly , extreme HPD widths are some what reduced in these cases relati ve to Cases Ia and IV a. T aken together , the HPD set results conﬁrm that redshift uncertainty is sensiti ve to the sig- nal intensity prior, though the magnitude of this sensitivity is less dramatic than its effect on detection decisions and the resulting marginal MAP estimates for redshift. S.7.4. Redshift hyperparameter sensitivity . F or redshift, we consider four scaled Beta priors ov er the interv al (0 , 4] : • Case Ib - Science-informed: a 3 = 3 , b 3 = 3 • Case IIb - Uniform: a 3 = 1 , b 3 = 1 B A YESIAN EMISSION LINE DETECTION 43 Case IIIb Case IVb Case Ib Case IIb 0 1 2 3 4 0 1 2 3 4 0 10 20 0 10 20 0 10 20 0 10 20 30 Redshift Count F I G S . 1 5 . Each histogram depicts the distribution of mar ginal MAP redshift estimates for a differ ent hyperpa- rameter setting for r edshift parameter . 0.001 0.010 0.100 1.000 Case Ib Case IIb Case IIIb Case IVb Hyper parameter combination T otal HPD width F I G S . 1 6 . Each histogr am depicts the total HPD width (on the log scale) for a differ ent hyperparameter setting for r edshift parameter . • Case IIIb - T ightly concentrated symmetric: a 3 = 10 , b 3 = 10 • Case IVb - Ske wed: a 3 = 3 , b 3 = 1 . S.7.4.1. Impact on log posterior odds and FPRs. Figure S.14 shows histograms of the log posterior odds under the four prior speciﬁcations. The overall shapes of the distributions 44 are qualitati vely similar , though the number of spectra identiﬁed as containing emission lines v aries from 203 to 208 across settings. T able S.3 reports the rejection rates on JWST data and corresponding FPRs on simulated data. Both quantities e xhibit some v ariation across prior choices, but the magnitude of sensi- ti vity is small relative to that observed for the background or signal intensity hyperparame- ters. In particular , the largest FPR across the four cases remains below 0.047, indicating that emission line detection is relati vely stable under these prior v ariations. S.7.4.2. Impact on r edshift distributions. Figure S.15 shows histograms of the marginal MAP redshift estimates across the 213 spectra. The primary differences across prior settings occur in the tails of the redshift estimate distribution, particularly near zero and four . Case IIIb places essentially no prior mass near the boundaries, and accordingly yields no MAP es- timates near zero or four . The remaining cases permit more e xtreme redshift values. Notably , the distribution under Case IIIb resembles the empirical distribution of astronomer-reported redshifts more closely than the other settings. Overall, aside from the boundaries, the distri- bution of the majority of g alaxies remains broadly similar across prior speciﬁcations. S.7.4.3. Impact on redshift HPD set widths. Figure S.16 displays boxplots of total 99.865% HPD set widths (on the log scale) for redshift under each prior . Sensitivity of HPD widths to the redshift prior is minimal. The boxplots are nearly indistinguishable across cases, with dif ferences appearing primarily in values near the boundaries. Cases Ib, IIb, and IVb each contain at least one spectrum with total HPD width exceeding one, while Case IIIb ex- hibits slightly smaller extremes, consistent with its reduced prior mass near the boundaries (that is, there are fe wer plausible redshift values under this prior). Overall, the HPD set widths are largely stable to these changes in prior h yperparameters.

A scalable Bayesian framework for galaxy emission line detection and redshift estimation

Original Paper

Comments & Academic Discussion

Leave a Comment

Original Paper

Related Papers

Comments & Academic Discussion

Leave a Comment