astro-ph.IMJun 13, 2025

Statistical Machine Learning for Astronomy -- A Textbook

Authors: Yuan-Sen Ting

Organizations: The Ohio State University

Read PDF
Ask AI
Cite

Copy citation

DOINot available for this paper

Abstract

This textbook gives a systematic treatment of statistical machine learning for astronomical research, based on Bayesian inference. It shows how modern data analysis techniques connect to classical statistical methods and follow from familiar statistical principles. The Bayesian perspective, used throughout, gives priority to uncertainty quantification and statistical rigor, which scientific inference in astronomy needs. The textbook progresses from probability theory and Bayesian inference through supervised learning, including linear regression with measurement uncertainties, logistic regression and classification. Unsupervised learning topics cover Principal Component Analysis and clustering methods. We then introduce sampling and Markov Chain Monte Carlo, followed by Gaussian Processes as probabilistic nonparametric methods and neural networks in the broader statistical context. Our theory-focused approach derives each method from first principles with complete mathematical development and then shows its astronomical applications. We put statistical insight first and explain why algorithms work, when they are appropriate, and how they connect to broader statistical principles. From a solid grounding in classical methods and their theory, the treatment leads to modern techniques such as neural networks. With this grounding, readers can apply the methods to astronomical research with proper attention to assumptions, limits and uncertainty propagation, as large astronomical surveys require. The 22 companion tutorials apply the methods to astronomical data in Jupyter notebooks.

Figures & tables

Figure 2.1: Illustration of deterministic versus random variables. The dashed vertical line represents a true (deterministic) value of some quantity. The blue curve shows the probability distribution of measured values when accounting for uncertainties. Individual red dots represent specific measurements, which scatter around the true value due to various sources of uncertainty. This demonstrates how a fundamentally deterministic quantity becomes a random variable when we consider measurement or other uncertainties.
Figure 2.2: Visualization of a standard normal (Gaussian) distribution with mean μ=0\mu=0 and standard deviation σ=1\sigma=1 . The dashed vertical line indicates the mean ( μ\mu ), which corresponds to the peak of the distribution. The green arrow shows one standard deviation ( σ\sigma ) from the mean, a key measure of the distribution’s spread. The Gaussian distribution is symmetric around its mean, with approximately 68% of the probability mass falling within one standard deviation of the mean. This distribution naturally emerges when many small, independent effects combine additively, making it ubiquitous in astronomical measurements due to the Central Limit Theorem.
Figure 2.3: Visualization of a Poisson distribution with mean rate parameter λ=5\lambda=5 . The discrete probability mass function is shown as blue bars, with the mean value indicated by the dashed vertical line. A key property of the Poisson distribution is that its standard deviation equals λ\sqrt{\lambda} (shown by the orange arrow), making the distribution relatively more spread out for smaller values of λ\lambda and more concentrated for larger values. This distribution naturally describes counting processes in astronomy, such as photon arrivals, cosmic ray hits, or galaxy counts in small volumes, where events occur independently at a constant average rate.
Figure 2.4: Illustration of how the Poisson distribution (blue bars) approaches a Gaussian distribution (red line) as the mean rate parameter λ\lambda increases. For small λ\lambda (e.g., λ=1\lambda=1 ), the distribution is highly discrete and skewed. As λ\lambda increases ( λ=5,10\lambda=5,10 ), the distribution becomes more symmetric and the discrete nature is less pronounced. By λ=50\lambda=50 , the Poisson distribution is nearly indistinguishable from a Gaussian distribution with the same mean and variance ( σ2=λ\sigma^{2}=\lambda ). This convergence, a consequence of the Central Limit Theorem, explains why photon counting statistics can often be approximated by Gaussian distributions when count rates are high.
Figure 2.5: Visualization of a power law distribution using the Salpeter Initial Mass Function (IMF) as an example. The Initial Mass Function (IMF) describes the distribution of stellar masses at birth—essentially how many stars of each mass are formed in a given population. The distribution follows ξ(M)∝M−α\xi(M)\propto M^{-\alpha} with α=2.35\alpha=2.35 (Salpeter’s value), plotted on logarithmic scales. The vertical dashed lines indicate the physical limits of stellar masses ( MminM_{\mathrm{min}} and MmaxM_{\mathrm{max}} ). This power law behavior emerges from the scale-invariant nature of the star formation process, where the same physical mechanisms operate across many orders of magnitude in mass. The negative slope indicates that lower-mass stars are much more numerous than high-mass stars, with the frequency decreasing by a factor of ∼\sim 102.35=22210^{2.35}=222 for every factor of 10 increase in mass. The distribution is normalized over the finite mass range, though it appears as a simple straight line in log-log space.
Figure 2.6: Illustration of joint and conditional probability distributions using stellar temperature and luminosity. The left panel shows the joint distribution p(T,L)p(T,L) describing the probability density of finding stars with different combinations of temperature and luminosity. The red plane indicates a slice through the joint distribution at a fixed temperature T=6000T=6000 K. The right panel shows the corresponding conditional distribution p(L∣T=6000K)p(L|T=6000\text{K}) , representing the probability distribution of stellar luminosities for stars with temperature T=6000T=6000 K.
Figure 2.7: Visualization of marginalization for a bivariate Gaussian distribution. The central panel shows contours of a joint probability distribution p(x,y)p(x,y) with positive correlation ( ρ=0.7\rho=0.7 ). The top and right panels show the marginal distributions p(x)p(x) and p(y)p(y) respectively, obtained by integrating the joint distribution over the other variable. The green arrows indicate the marginalization process: integrating over yy (vertical arrow) gives the marginal distribution of xx (top panel), while integrating over xx (horizontal arrow) gives the marginal distribution of yy (right panel). This illustrates how marginalization reduces a multivariate distribution to its component distributions by accounting for all possible values of the marginalized variables.
Figure 2.8: Visualization of dependent versus independent variables. Left panels show a case where the variables xx and yy are independent — note how the conditional distributions p(y∣x)p(y|x) remain identical regardless of the value of xx . This means knowing xx provides no information about yy . Right panels show dependent variables where p(y∣x)p(y|x) changes with different values of xx , demonstrating that knowledge of xx informs our expectations about yy . Top panels show the conditional distributions for two different values of xx (solid blue line: x=1x=1 , dashed red line: x=−1x=-1 ). Bottom panels show the corresponding joint distributions with vertical lines indicating where the conditional slices were taken.
OutcomeFairBadWorseWorst
61/61/81/121/20
51/61/81/121/20
41/61/61/81/10
31/61/61/61/8
21/65/241/41/5
11/65/241/31/2
Table 2.1: Probability distributions for different die models. Each column represents a different hypothesis about the die’s behavior, with Fair representing an unbiased die and successive models showing increasing bias toward lower numbers.
Figure 2.9: Top panel: Likelihood functions p(x∣μ)p(x|\mu) for different proposed values of μ\mu . When μ=2\mu=2 (red curve), the probability of measuring x=6x=6 is very low. When μ=5\mu=5 (green curve), this probability is much higher. Bottom panel: The likelihood p({6}∣μ)p(\{6\}|\mu) as a function of μ\mu , showing how likely each proposed mean value is given our single observation. This likelihood is proportional to our posterior distribution (assuming an uninformative prior). The maximum occurs at μ=6\mu=6 , which makes intuitive sense—with a single measurement and no informative prior, our best estimate for the mean would be exactly what we measured.
Figure 2.10: Top panel: Individual likelihood functions for different values of μ\mu , with vertical lines marking our two observations at x=4x=4 and x=6x=6 . Bottom panel: The joint likelihood—showing how the probability of observing both measurements varies with μ\mu . The maximum likelihood estimate occurs at μ=5\mu=5 , the average of our two measurements.
Figure 3.1: Visualization of how different moments weight a probability distribution. Top panel: The black curve shows a normalized probability distribution p(x)p(x) , while the colored curves show xkp(x)x^{k}p(x) for different values of kk . As kk increases, the peaks of xkp(x)x^{k}p(x) shift toward larger values of xx , demonstrating how higher moments become increasingly sensitive to the tails of the distribution. Bottom panel: The relative contribution of different powers of xx ( xkx^{k} ) to each moment, showing how higher moments ( k>1k>1 ) give progressively more weight to larger values. This illustrates why higher moments are particularly sensitive to extreme values in a distribution—a feature especially relevant in astronomy where rare, extreme objects often carry important physical information.
Figure 3.2: Visualization of how probability distributions transform under linear transformations. Top panel: The blue curve shows the original probability distribution pX(x)p_{X}(x) and the green curve shows the transformed distribution pY(y)p_{Y}(y) under the linear transformation Y=2X+1Y=2X+1 . The shaded regions have equal areas, demonstrating probability conservation. Bottom panel: Illustration of how a uniform grid transforms under the same linear transformation. The spacing between grid lines increases by a factor of ∣a∣=2|a|=2 , necessitating a corresponding decrease in probability density by a factor of 1/∣a∣1/|a| to preserve total probability. This visualization demonstrates why the Jacobian factor ∣dx/dy∣|dx/dy| appears in probability transformation formulas.
Figure 3.3: Visualization of expectation’s linearity property using a bivariate Gaussian distribution. Left panel: Original distribution with expectation E[X]\mathbb{E}[X] marked by the red cross. Right panel: The transformed distribution under Y=2X+1Y=2X+1 , showing how the expectation transforms linearly: E[Y]=2E[X]+1\mathbb{E}[Y]=2\mathbb{E}[X]+1 . The contours represent probability density levels, demonstrating how the entire distribution stretches and shifts while maintaining its shape. The red dashed lines highlight how the center (expectation) transforms according to the same linear rule, regardless of the distribution’s specific shape or spread.
Figure 3.4: Visualization of expectation’s additivity property using bivariate Gaussian distributions. Left panel: Original distributions of random variables XX (blue contours) and YY (green contours) with their respective expectations E[X]\mathbb{E}[X] and E[Y]\mathbb{E}[Y] marked by colored crosses. The contours represent probability density levels of each distribution, and the dashed lines mark the coordinate positions of the expectations. Right panel: The distribution of the sum Z=X+YZ=X+Y (purple contours) with its expectation E[Z]=E[X+Y]\mathbb{E}[Z]=\mathbb{E}[X+Y] marked by the red cross. The figure demonstrates that E[X+Y]=E[X]+E[Y]\mathbb{E}[X+Y]=\mathbb{E}[X]+\mathbb{E}[Y] , showing how expectations add component-wise regardless of the correlation structure in the original distributions.
Figure 3.5: Visualization of the law of total expectation using stellar mass-luminosity relationships. Left panel: Joint distribution of stellar mass and luminosity (in solar units and log scale), shown as blue contours. The red cross marks the overall expectation E[M]\mathbb{E}[M] , while colored circles show conditional expectations E[M∣L]\mathbb{E}[M|L] at three different luminosities (marked by dashed lines). Right panel: The conditional expectation E[M∣L]\mathbb{E}[M|L] as a function of luminosity (blue line), demonstrating how we can estimate stellar masses given their luminosities. The gray shaded region shows the marginal distribution of luminosity p(L)p(L) , which accounts for observational selection effects in magnitude-limited surveys. The red dashed line marks the overall expectation E[M]\mathbb{E}[M] , illustrating how integrating E[M∣L]\mathbb{E}[M|L] weighted by p(L)p(L) recovers the average stellar mass E[M]\mathbb{E}[M] , even when direct mass measurements are challenging.
Figure 3.6: Visualization of how variance and standard deviation transform under linear transformations. Left panel: Original bivariate distribution (blue contours) with its expectation marked by the red cross and variance structure indicated by the red dashed ellipse. The ellipse represents a constant probability contour of the distribution, with its shape and orientation determined by the covariance matrix. Right panel: The transformed distribution (green contours) under Y=2X+1Y=2X+1 . While the mean shifts by the same linear transformation ( E[Y]=2E[X]+1\mathbb{E}[Y]=2\mathbb{E}[X]+1 ), the standard deviation scales linearly ( σY=2σX\sigma_{Y}=2\sigma_{X} ) while the variance scales quadratically ( Var[Y]=4Var[X]\text{Var}[Y]=4\text{Var}[X] ), as shown by the enlarged red dashed ellipse.
Figure 3.7: Visualization of how correlation shapes bivariate normal distributions. The panels show three distributions with identical marginal variances but different correlation structures. Left panel: Uncorrelated variables ( ρ=0\rho=0 ) produce circular contours, with the principal axes of the red dashed ellipse aligned with the coordinate axes. Middle panel: Positive correlation ( ρ=0.7\rho=0.7 ) stretches the distribution along the diagonal (red dashed line), indicating that higher values of one variable tend to occur with higher values of the other. Right panel: Negative correlation ( ρ=−0.7\rho=-0.7 ) tilts the distribution in the opposite direction, showing that high values of one variable tend to occur with low values of the other. The red dashed ellipses, representing contours of constant probability, rotate and deform based on the correlation strength, while maintaining the same total variance (area). This illustrates how correlation captures the directionality of the relationship between variables without affecting their individual scales.
Figure 3.8: Visualization of how variances add for independent random variables. Left panel: Distribution of XX (blue contours) with different variances in each dimension (1.0, 0.5). Middle panel: Distribution of YY (green contours) with different variances (0.5, 1.0). Right panel: Distribution of their sum X+YX+Y (purple contours). The red dashed ellipses show the variance structure of each distribution. For independent variables, the variances add component-wise: the variance of the sum in each dimension equals the sum of the individual variances in that dimension. This is evident in the right panel, where Var[X+Y]=(1.5,1.5)\text{Var}[X+Y]=(1.5,1.5) results from adding the corresponding variances. This increase in variance is intuitive: when we add two independent uncertain quantities, our uncertainty naturally grows since each variable contributes its own source of randomness to the final sum.
Figure 3.9: Demonstration of how sample mean uncertainty decreases with sample size. Left panel: The black curve shows the original distribution of individual measurements (with true mean shown by red dashed line). The colored histograms show the distributions of sample means for different sample sizes ( N=5N=5 , 20, and 100). As NN increases, the distribution of sample means becomes increasingly concentrated around the true mean, demonstrating how larger samples provide more precise estimates. Right panel: Quantitative analysis of how the standard error (uncertainty in the sample mean) scales with sample size. The red dashed line shows the theoretical prediction σ/N\sigma/\sqrt{N} , while blue points show the measured uncertainties from numerical simulations.
Figure 3.10: Illustration of the bootstrap resampling process and its ability to estimate sampling variance of the first moment. Left panel: Original sample {x1,...,xN}\{x_{1},...,x_{N}\} (black tick marks) drawn from true distribution p(x)p(x) (black dashed curve), with three bootstrap resamples {x1∗,...,xN∗}\{x_{1}^{*},...,x_{N}^{*}\} (colored histograms) generated by sampling with replacement. Right panel: Distribution of bootstrap estimates μ^1∗=1N∑i=1Nxi∗\hat{\mu}_{1}^{*}=\frac{1}{N}\sum_{i=1}^{N}x_{i}^{*} from 1000 resamples (blue histogram) compared to theoretical sampling distribution of first moment estimator (red dashed curve). The agreement between bootstrap variance Var(μ^1∗)\text{Var}(\hat{\mu}_{1}^{*}) and theoretical variance σ2/N\sigma^{2}/N demonstrates bootstrap’s accuracy in estimating sampling variance without knowledge of p(x)p(x) . Black vertical line shows original sample estimate μ^1=1N∑i=1Nxi\hat{\mu}_{1}=\frac{1}{N}\sum_{i=1}^{N}x_{i} .
Figure 3.11: Visualization of the law of total variance using the stellar mass-luminosity relationship. The blue points show individual stars, demonstrating the scatter in stellar masses at each luminosity. The red line traces the conditional expectation EM[M∣L]\mathbb{E}_{M}[M|L] as a function of luminosity, with the red shaded region showing ±1σ\pm 1\sigma bounds of the conditional variance. The blue shaded region at L=0.5L⊙L=0.5L_{\odot} highlights EL[VarM[M∣L]]\mathbb{E}_{L}[\text{Var}_{M}[M|L]] , representing the average scatter in mass at fixed luminosity. The green shaded region and dashed lines on the right show the range of conditional expectations, whose variance VarL[EM[M∣L]]\text{Var}_{L}[\mathbb{E}_{M}[M|L]] captures how the average mass changes with luminosity. The law of total variance states that the total variance in stellar mass equals the sum of these two contributions: the average scatter at fixed luminosity ( EL[VarM[M∣L]]\mathbb{E}_{L}[\text{Var}_{M}[M|L]] ) plus the variance in the mass-luminosity relationship itself ( VarL[EM[M∣L]]\text{Var}_{L}[\mathbb{E}_{M}[M|L]] ). This decomposition helps us understand how much of the mass variance comes from intrinsic scatter versus systematic trends with luminosity.
Figure 4.1: Visualization of the likelihood function for linear regression with Gaussian noise. The cyan curves show model predictions y(x,w)y(\mathbf{x},\mathbf{w}) for two different parameter values w\mathbf{w} . Yellow distributions indicate the Gaussian uncertainty in possible observations at each prediction point, while red dots show actual observations. (a) A model with appropriate linear coefficients, where observations scatter naturally within the expected uncertainty. (b) A model with incorrect coefficients, showing the deviations of the model predictions from the data that cannot be explained by data measurement uncertainty alone.
Figure 4.2: Demonstration of how different choices of uncertainty values affect the likelihood in maximum likelihood estimation. Panel (a): With underestimated uncertainties ( σ=0.1\sigma=0.1 ), the narrow Gaussian distributions mean data points far from the model contribute very small likelihood terms, resulting in a low total likelihood. Panel (b): When uncertainties are correctly estimated ( σ=0.3\sigma=0.3 ), we achieve the maximum likelihood—this represents the optimal balance between the width of the Gaussian distributions and their heights. Panel (c): With overestimated uncertainties ( σ=0.6\sigma=0.6 ), while the Gaussian distributions are wide enough to encompass most points, the height of each Gaussian (which scales as 1/σ1/\sigma ) becomes smaller, resulting in smaller likelihood values overall.
Figure 4.3: Schematic illustration showing the globally convex log-likelihood surface in (w,β)(w,\beta) space. The bowl-shaped contours demonstrate why optimization can be achieved in exactly two steps: since the optimal ww depends only on relative residual weights (independent of β\beta scale), we can first optimize ww for any initial β\beta , then directly find the optimal β\beta using that ww . This two-step convergence is guaranteed by the convex geometry, regardless of starting point.
Figure 4.4: Demonstration of polynomial fitting with increasing model complexity. Each panel shows the true underlying function (blue dashed line), training data points (red), test data points (green), and the fitted polynomial (pink) of different orders. As the polynomial order increases from 1 to 15, we observe how the model’s flexibility changes. The linear fit (order 1) is too rigid to capture the data’s pattern, while the high-order polynomial (order 15) follows the noise in the training data too closely. Intermediate orders provide a better balance between fitting the data and capturing the true underlying pattern.
Figure 4.5: Training and test error curves (on log scale) illustrating the three fundamental regimes in model complexity. The underfitting regime (red shading) shows high error on both training and test sets due to insufficient model flexibility. The good fit regime (green shading) represents optimal model complexity where both errors are minimized. The overfitting regime (blue shading) is characterized by very low training error but increasing test error as the model begins to fit noise in the training data. This behavior demonstrates the crucial balance between model complexity and generalization performance.
Figure 5.1: Illustration of why we need to go beyond simple maximum likelihood estimation. Both panels show the same underlying relationship (dashed blue line) fitted to observed data (orange points), but with different measurement uncertainties (shown as yellow distributions). Panel (a) shows a case with small measurement errors, where bootstrapped uncertainty estimates are reliable. Panel (b) shows the same fit applied to data with large measurement errors, where bootstrapping may provide misleading uncertainty estimates due to the limited sample size and complex error structure.
Figure 5.2: Visualization of Poisson-Gamma conjugacy in star count analysis. Panel (a) demonstrates how multiplying the Poisson likelihood (dotted blue line) with a Gamma prior (dashed green line) yields a Gamma posterior (solid red line, with shading). The prior Gamma( α=10\alpha=10 , β=2\beta=2 ) represents previous observations suggesting about 5 stars per region, which is updated by a new observation of 8 stars. Panel (b) shows how different prior parameters affect the posterior distribution, while keeping the same likelihood (8 stars observed). When the prior has higher precision ( α=20\alpha=20 , β=4\beta=4 , darker lines), the posterior stays closer to the prior mean. When the prior has lower precision ( α=5\alpha=5 , β=1\beta=1 , lighter lines), the posterior is more strongly influenced by the likelihood. This illustrates how the relative strength of prior knowledge versus new data shapes our updated beliefs.
Figure 5.3: Visualization of Gaussian conjugacy. Panel (a) demonstrates how multiplying a Gaussian prior (dashed green line) with a Gaussian likelihood (dotted blue line) yields a Gaussian posterior (solid red line, with shading). Panel (b) shows how different prior precisions affect the posterior distribution, while keeping the same likelihood. When the prior has high precision ( τ02=0.5\tau_{0}^{2}=0.5 , more concentrated), the posterior stays closer to the prior. When the prior has low precision ( τ02=4\tau_{0}^{2}=4 , more spread out), the posterior is dominated by the likelihood. This illustrates how the relative precisions of the prior and likelihood determine their influence on the posterior beliefs.
Figure 5.4: Visualization of the product of two 2D Gaussian distributions. The blue and orange dashed contours show the level sets of two original Gaussian distributions, representing our prior beliefs and data-driven information about linear regression parameters. Their means are marked by corresponding colored ‘+’ symbols. The filled red contours show the resulting posterior distribution obtained from their product. The red ‘+’ marks the posterior mean, which is a precision-weighted average of the input means. Note how the posterior distribution is more concentrated than either input distribution, reflecting the combination of information from both sources. The position and shape of the posterior captures how our prior knowledge and observed data combine to yield updated beliefs.
Figure 5.5: Comparison of maximum likelihood and Bayesian approaches to linear regression. Panel (a) shows the maximum likelihood estimate (MLE), which provides a single best-fit line through the data. The error bars on the data points indicate the assumed measurement uncertainty ( σ\sigma ). Panel (b) shows the Bayesian approach with the same data, where the light blue lines represent different possible fits drawn from the posterior distribution over regression parameters. Each line represents a different plausible fit given our data and prior assumptions. The solid blue line shows the posterior mean prediction, and the shaded region represents the 95% credible interval. Both approaches identify the same underlying trend, but the Bayesian approach naturally quantifies our uncertainty in the fit, with credible intervals that widen away from the data points where we have less constraint on the model parameters.
Figure 5.6: Effect of regularization strength ( λ=σ2/η2\lambda=\sigma^{2}/\eta^{2} ) in Bayesian linear regression. Panel (a) shows weak regularization (small λ\lambda , corresponding to precise measurements or uncertain prior), where the model has more freedom to fit the data, resulting in wider prediction bands that reflect greater parameter uncertainty, especially away from the data points. Panel (b) shows strong regularization (large λ\lambda , corresponding to noisy measurements or confident prior), where the prior more strongly constrains the parameters toward smaller values, leading to narrower prediction bands and simpler, more conservative fits. The blue lines show samples from the posterior distribution of possible fits, while the shaded region indicates the 95% credible interval for predictions. This demonstrates how the regularization parameter λ\lambda controls the trade-off between model flexibility and constraint—a key mechanism for preventing overfitting.
Figure 5.7: Visualization of the posterior predictive distribution in Bayesian linear regression. The red points show observed data with error bars indicating measurement uncertainties. The light blue lines represent individual regression fits sampled from the posterior distribution p(w∣D)p(\mathbf{w}|\mathcal{D}) , illustrating the range of models consistent with our data and prior assumptions. The shaded blue region shows the 95% credible interval for predictions. At three specific input points x∗x_{*} , we visualize the complete predictive distribution p(y∗∣x∗,D)p(y_{*}|x_{*},\mathcal{D}) as probability density curves. The mean of each predictive distribution corresponds to our best prediction, while the width reflects our total uncertainty, which increases away from the data points where we have less constraint on the model parameters.
Figure 6.1: Illustration of how measurement uncertainties affect regression. Both panels show the same underlying relationship (dashed blue line) with true positions (hollow circles) and their observed values (solid orange points) connected by gray arrows. Panel (a) shows the case with uncertainties only in the dependent variable ( yy ), where the vertical scatter of observed points around their true positions does not systematically bias our inference. Panel (b) shows what happens when both variables have uncertainties: observed points now scatter in both directions from their true positions, artificially inflating the apparent range of our independent variable ( xx ). This horizontal scatter systematically weakens any relationship we try to measure, leading to attenuation bias. The error bars in both panels represent the measurement uncertainties ( σy\sigma_{y} in panel a, σx\sigma_{x} and σy\sigma_{y} in panel b), illustrating how the observed points deviate from their true positions due to measurement noise.
Figure 6.2: Demonstration of how measurement uncertainties affect the apparent dynamic range of our data. Both panels show the same underlying relationship (solid black line) and true data points (hollow green circles) scattered by measurement uncertainties to their observed positions (solid orange circles). The green and orange shaded regions show the true and observed dynamic ranges respectively. Panel (a) shows the case where measurement uncertainty is small compared to the intrinsic spread in x-values ( σx/σrange=0.1\sigma_{x}/\sigma_{\text{range}}=0.1 ), resulting in a modest increase in the observed range. Panel (b) demonstrates the case where measurement uncertainty equals the intrinsic spread ( σx/σrange=1.0\sigma_{x}/\sigma_{\text{range}}=1.0 ), causing an inflation of the observed dynamic range. This artificial expansion of the x-range without a corresponding increase in the y-range directly leads to attenuation bias in our parameter estimates, as the observed distribution appears more stretched horizontally than it truly is.
Figure 6.3: Demonstration of attenuation bias in parameter estimation. Both panels show the same underlying relationship (solid black line) with observed data points (orange circles) and their fitted relationships (dashed green lines). Panel (a) shows the case of small measurement uncertainty ( σx/σrange=0.1\sigma_{x}/\sigma_{\text{range}}=0.1 ), where the fitted slope is only slightly attenuated ( λ≈0.99\lambda\approx 0.99 ). Panel (b) shows the case of large measurement uncertainty ( σx/σrange=1.0\sigma_{x}/\sigma_{\text{range}}=1.0 ), where the attenuation becomes severe ( λ≈0.50\lambda\approx 0.50 ). The attenuation factor λ\lambda , shown in the legend of each panel, quantifies how much the observed slope underestimates the true slope. This demonstrates that when measurement uncertainties become comparable to the intrinsic spread in our independent variable, our ability to recover the true relationship is compromised.
Figure 6.4: Demonstration of how measurement uncertainties affect parameter estimation through attenuation bias. The main plot shows how the attenuation factor λ\lambda depends on the ratio of measurement uncertainty to intrinsic spread ( σx/σrange\sigma_{x}/\sigma_{\text{range}} ). The solid black line shows the theoretical prediction λ=1/(1+(σx/σrange)2)\lambda=1/(1+(\sigma_{x}/\sigma_{\text{range}})^{2}) . Three inset panels demonstrate the effect for different uncertainty ratios: small ( σx/σrange=0.1\sigma_{x}/\sigma_{\text{range}}=0.1 ), moderate ( σx/σrange=0.5\sigma_{x}/\sigma_{\text{range}}=0.5 ), and large ( σx/σrange=1.0\sigma_{x}/\sigma_{\text{range}}=1.0 ). In each inset, the solid black line shows the true relationship, orange points show observed data (with both x and y uncertainties), and the dashed green line shows the attenuated relationship recovered from the data. This illustrates how measurement uncertainties comparable to or larger than the intrinsic spread lead to systematically shallower slopes in our parameter estimates.
Figure 6.5: Demonstration of how we estimate true positions in a hierarchical model. Both panels show the same underlying relationship (dashed blue line) and observed data points with measurement uncertainties (error bars). Panel (a) shows the raw observed data, where both x and y measurements carry uncertainties. Panel (b) visualizes our inference process: for each observation, we estimate both its true x-position (blue points) and corresponding model-predicted y-value. Gray arrows show the two-step correction: first a horizontal shift from the observed x to its estimated true value, then a vertical shift to the model prediction. The magnitude of these corrections depends on the measurement uncertainties ( σx\sigma_{x} and σy\sigma_{y} ) and the model parameter ( β\beta ). The estimated true positions (blue points) lie exactly on the model line because they represent our complete model prediction (xtrue,βxtrue)(\mathbf{x}_{\text{true}},\beta\mathbf{x}_{\text{true}}) , while the intermediate points (green) show the state after x-correction but before y-correction. This hierarchical estimation process, updating both xtrue\mathbf{x}_{\text{true}} values and β\beta , forms the basis of Deming regression.
Figure 6.6: Demonstration of how Deming regression corrects for attenuation bias compared to ordinary least squares (OLS, also known as the maximum likelihood solution when there is no xx -uncertainty). Panel (a) shows simulated data (orange points with error bars) from a true linear relationship y=3xy=3x (solid black line) with measurement uncertainties σx=1.0\sigma_{x}=1.0 and σy=0.5\sigma_{y}=0.5 . The OLS fit (green dashed line) systematically underestimates the true slope due to attenuation bias, while Deming regression (blue dashed line) better recovers the true relationship by properly accounting for measurement uncertainties in both variables. Panel (b) demonstrates how these estimators behave as the relative measurement uncertainty ( σx/σy\sigma_{x}/\sigma_{y} ) varies: OLS estimates become increasingly attenuated as xx -uncertainty grows, while Deming regression maintains better accuracy by appropriately weighting the uncertainties. When xx -uncertainty is small compared to yy -uncertainty ( σx≪σy\sigma_{x}\ll\sigma_{y} ), both methods converge to the true slope as standard regression assumptions become valid.
Figure 6.7: Demonstration of heteroscedastic Deming regression and its numerical solution. Panel (a) shows simulated data (orange points) with heterogeneous measurement uncertainties (gray error bars) drawn from uniform distributions: σx∼U(0.5,1.5)\sigma_{x}\sim U(0.5,1.5) and σy∼U(0.3,0.7)\sigma_{y}\sim U(0.3,0.7) . The true relationship y=3xy=3x (solid black line) is compared with ordinary least squares (green dashed line) and heteroscedastic Deming regression (blue dashed line). Panel (b) illustrates the numerical solution method: the objective function F(β)F(\beta) (solid blue line) and its absolute value ∣F(β)∣|F(\beta)| (dashed red line) are plotted against candidate slope values. The Deming regression slope is found where F(β)F(\beta) crosses zero, or equivalently where ∣F(β)∣|F(\beta)| reaches its minimum.
Figure 7.1: Demonstration of binary classification in a two-dimensional feature space. Stars (blue) and galaxies (red) form distinct clusters based on two observable features. The dashed green line shows the decision boundary learned by logistic regression, which separates the two classes. A new object (yellow star) is classified based on its position relative to this boundary, with the model assigning probabilities to each class membership. This illustrates how classification algorithms learn to partition the feature space and make predictions for new observations. The logistic regression—a linear classifier—assumes a linear decision boundary as shown. The features could represent various astronomical measurements that help distinguish stars from galaxies, such as photometric colors, morphological parameters (comparing PSF vs model magnitudes), or astrometric measurements from missions like Gaia (proper motions and parallaxes).
Figure 7.2: Visualization of generative versus discriminative approaches in classification. Panel (a) shows the generative model perspective: learning the class-conditional densities p(x∣Ck)p(x|C_{k}) for each class separately. This illustrates how generative models capture the full complexity of feature distributions within each class, enabling tasks beyond classification like anomaly detection (e.g., identifying quasars that don’t match either stellar or galactic distributions) and synthetic data generation. Panel (b) shows the discriminative model perspective: directly learning the posterior probabilities p(Ck∣x)p(C_{k}|x) of class membership. The vertical purple line indicates the decision boundary where the probabilities of both classes are equal ( p(C1∣x)=p(C0∣x)=0.5p(C_{1}|x)=p(C_{0}|x)=0.5 ). While simpler to learn, the discriminative approach focuses solely on classification, potentially missing valuable information about the underlying physical distributions that generative models preserve.
Figure 7.3: Illustration of the problem of simply treating classification as regression. Panel (a) shows the 2D feature space with two classes (cyan and yellow) and the linear decision boundary (purple dashed line) where y=0.5y=0.5 . The weight vector w\mathbf{w} (red arrow) is perpendicular to the decision boundary, determining its orientation. Gray lines show contours of constant yy values (-1, 0, 1, 2), demonstrating how the linear predictor continues to grow indefinitely as we move away from the boundary. Panel (b) shows the same scenario in 3D, where the vertical axis represents the raw output y=wTxy=\mathbf{w}^{T}\mathbf{x} of the linear model. While the decision boundary at y=0.5y=0.5 appears reasonable, the unbounded nature of yy means that points far from the boundary receive predictions well outside the valid probability range of [0,1], leading to potential issues in model fitting.
Figure 7.4: Illustration of how direct linear regression can lead to suboptimal decision boundaries in classification. Panel (a) shows the ideal decision boundary (purple line) that correctly separates two classes when considering only the main clusters, with gray lines showing contours of constant yy values from the linear predictor wTx\mathbf{w}^{T}\mathbf{x} . Panel (b) demonstrates how the boundary can be skewed (green line) when attempting direct linear regression with outliers: because the linear predictor grows unboundedly, outlying points produce very large errors ( y≫1y\gg 1 ) that dominate the least squares optimization. The model shifts the boundary to reduce these large errors, even though this compromises the classification of the main clusters. This illustrates why we need a tapering function like the sigmoid to ensure predictions remain bounded between 0 and 1, preventing outliers from having undue influence on the boundary placement.
Figure 7.5: The sigmoid function and its role in logistic regression. Panel (a) shows the sigmoid function σ(z)=1/(1+e−z)\sigma(z)=1/(1+e^{-z}) , which “squashes” any real number into the range [0,1]. When z≪0z\ll 0 , the function approaches 0; when z≫0z\gg 0 , it approaches 1; and at z=0z=0 , it takes the value 0.5. Panel (b) shows how this creates probability contours in feature space: the decision boundary (purple line) corresponds to σ(z)=0.5\sigma(z)=0.5 , while the gray contour lines show both the linear predictor zz values and their corresponding sigmoid probabilities σ(z)\sigma(z) . This demonstrates how the sigmoid function ensures predictions remain bounded between 0 and 1 regardless of distance from the decision boundary.
Figure 7.6: Visualization of how logistic regression transforms unbounded linear predictions into probabilities. Panel (a) shows the linear predictor log⁡(p(C1∣X)/p(C0∣X))=wTx+b\log(p(C_{1}|X)/p(C_{0}|X))=\mathbf{w}^{T}\mathbf{x}+b in feature space, where the surface extends infinitely in the vertical direction. The purple dashed line indicates the decision boundary where the log-odds ratio equals zero. Panel (b) shows the same predictor after applying the sigmoid function, transforming the unbounded log-odds into probabilities p(C1∣X)p(C_{1}|X) bounded between 0 and 1. The sigmoid function “squashes” the linear predictor while preserving the location of the decision boundary (where p(C1∣X)=0.5p(C_{1}|X)=0.5 ). Blue crosses and red circles represent two classes of data points, demonstrating how points far from the boundary receive probability predictions close to 0 or 1, while points near the boundary receive intermediate probabilities.
Figure 7.7: Illustration of how Gaussian generative models with equal covariance matrices naturally lead to linear decision boundaries. The contour lines show the probability density of two Gaussian distributions with means μ0\mu_{0} and μ1\mu_{1} (shown as larger markers) sharing the same covariance matrix Σ\Sigma . Sample points (crosses and circles) are drawn from each distribution. The optimal decision boundary (purple dashed line) emerges naturally from the equality of the covariance matrices and is perpendicular to the direction w∝Σ−1(μ1−μ0)\mathbf{w}\propto\Sigma^{-1}(\mu_{1}-\mu_{0}) . This demonstrates how our assumptions about the generative model’s covariance structure directly determine the form of the optimal classifier.
Figure 7.8: Illustration of how the covariance structure affects the optimal decision boundary in Gaussian generative models. Both panels share the same class means μ0\mu_{0} and μ1\mu_{1} but have different covariance matrices Σ\Sigma . Panel (a) shows the case of correlated features, where off-diagonal terms in Σ\Sigma create tilted elliptical contours, resulting in a decision boundary that accounts for this correlation. Panel (b) demonstrates uncorrelated features with different variances in each direction, producing axis-aligned elliptical contours of varying spread. In each case, the weight vector w∝Σ−1(μ1−μ0)\mathbf{w}\propto\Sigma^{-1}(\mu_{1}-\mu_{0}) (black arrows) points in different directions because Σ−1\Sigma^{-1} transforms the mean difference vector according to the feature covariance structure. The shape of the contours directly visualizes the role of Σ\Sigma in determining how the features vary, while the decision boundaries (purple dashed lines) show how Σ−1\Sigma^{-1} affects the optimal classification rule.
Figure 7.9: Demonstration of how class priors affect the decision boundary in logistic regression. Both panels show the same feature distributions (Gaussian contours) and weight vector orientation (determined by w=Σ−1(μ1−μ0)\mathbf{w}=\Sigma^{-1}(\mu_{1}-\mu_{0}) ), but with different class priors. Panel (a) shows the case of equal class priors ( ln⁡(P(C1)/P(C0))=0\ln(P(C_{1})/P(C_{0}))=0 ), where the decision boundary (purple dashed line) is placed solely based on the feature distributions. Panel (b) shows the effect when Class 1 is ten times more likely than Class 0 ( ln⁡(P(C1)/P(C0))=ln⁡(10)\ln(P(C_{1})/P(C_{0}))=\ln(10) ), causing the boundary to shift towards Class 0 to reflect this prior knowledge. The number of sample points (crosses for Class 0, circles for Class 1) also reflects this class imbalance. This illustrates how the prior term in the intercept bb can independently shift the decision boundary without affecting its orientation, allowing the classifier to account for known class imbalances in the training data.
Figure 7.10: Demonstration of how feature transformations can improve classification performance. Panel (a) shows the original data in Cartesian coordinates, where two classes (blue crosses for Class 0, red circles for Class 1) form curved, non-linearly separable patterns. The purple dashed line shows the best linear decision boundary, which poorly separates the classes. Panel (b) shows the same data transformed into polar coordinates ( θ\theta -r space), where the classes become linearly separable. The contours show the Gaussian distributions fit to each class in this transformed space, and the purple dashed line shows the optimal linear decision boundary. This illustrates how appropriate feature transformations can make logistic regression’s linear boundary assumption more valid by reshaping the data distribution.
Figure 7.11: Visualization of gradient descent optimization in weight space. The contours represent the loss landscape E(w)E(\mathbf{w}) , with darker colors indicating higher loss values. The red path shows the trajectory of gradient descent starting from a high-loss region (red star), with arrows indicating the negative gradient direction −∇E(w)-\nabla E(\mathbf{w}) at each step. The learning rate η\eta determines the step size along this path. The landscape illustrates common optimization challenges: multiple local minima that could trap the optimization, varying gradient magnitudes across the space, and plateau regions where gradients become small. The global minimum (annotated) represents the optimal weight configuration, while local minima demonstrate potential convergence points that may not be globally optimal.
Figure 7.12: Visualization of gradient contributions in logistic regression. The plot shows data points from two classes (blue circles for class 0, orange squares for class 1) along with their gradient vectors. The dashed line represents the current decision boundary. Each point’s gradient is represented by an arrow, with length proportional to both the prediction error and feature magnitude. Correctly classified points have shorter, grey arrows, while misclassified points have longer, colored arrows in their respective class colors. Note how the “Extreme” point, being far from the decision boundary, contributes a large gradient due to both its high prediction error and large feature magnitude. This visualization demonstrates two key properties of logistic regression gradients: (1) misclassified points contribute larger gradients than correctly classified ones, and (2) points with larger feature magnitudes generate proportionally larger gradients, as evidenced by the arrow lengths.
Figure 7.13: Comparison of Gradient Descent (GD) and Stochastic Gradient Descent trajectories in weight space. The contour plot shows the loss landscape E(w)E(\mathbf{w}) , with darker colors indicating lower loss values. The red path shows the deterministic GD trajectory, which follows the exact gradient at each step. The blue paths demonstrate multiple possible SGD trajectories, illustrating its stochastic nature. While GD follows a direct path toward a local minimum, SGD’s added noise allows it to explore more of the parameter space and potentially escape local minima. Starting from the same point (red star), GD converges to a local minimum (green star), while SGD paths show more exploration and can potentially find the global minimum.
Figure 7.14: Comparison of Stochastic Gradient Descent behavior with different learning rates η\eta . Both panels show the same loss landscape E(w)E(\mathbf{w}) with contours indicating level sets of the loss function, where darker colors represent lower loss values. Left: With a small learning rate ( η=0.05\eta=0.05 ), the trajectories show more conservative exploration, following the gradient more closely but potentially getting trapped in local minima. Right: With a larger learning rate ( η=0.2\eta=0.2 ), the paths exhibit more extensive exploration of the parameter space, potentially escaping local minima but risking overshooting optimal points. Multiple SGD paths (blue) starting from the same initialization point (red star) demonstrate how the combination of stochastic gradients and learning rate affects the optimization dynamics. In both cases, we observe the characteristic “noisy” trajectories of SGD compared to standard gradient descent, arising from the stochastic gradient estimates. This visualization highlights the crucial role of learning rate selection in balancing between exploitation and exploration.
Figure 7.15: Comparison of two learning rate scheduling strategies. Left: Inverse time decay schedule with η(t)=η01+γt\eta^{(t)}=\frac{\eta_{0}}{1+\gamma t} , showing a monotonic decrease in learning rate that starts rapid and becomes increasingly gradual. This schedule embodies the intuition that large steps are most beneficial early in training, with progressively finer adjustments needed as we approach convergence. Right: Cosine annealing schedule with restarts, where η(t)=ηmin⁡+12(η0−ηmin⁡)(1+cos⁡(tπ/T))\eta^{(t)}=\eta_{\min}+\frac{1}{2}(\eta_{0}-\eta_{\min})(1+\cos(t\pi/T)) . The cyclic nature of this schedule allows periodic returns to higher learning rates, potentially helping escape local minima. Each cycle begins with a phase of exploration (high learning rate) and ends with exploitation (low learning rate), similar to the repeated heating and cooling cycles in metallurgical annealing processes.
Figure 8.1: Illustration of pathological cases in One-vs-Rest (OvR) multi-class classification. The plot shows two binary decision boundaries: a solid red line separating class C1C_{1} from “not C1C_{1} ” and a dashed blue line separating class C2C_{2} from “not C2C_{2} ”. These boundaries create four distinct regions: a green region where both classifiers claim the point (both predict their positive class), a purple region where neither classifier claims the point (both predict their negative class, implying C3C_{3} by default), a red region where only C1C_{1} is predicted, and a blue region where only C2C_{2} is predicted. The ambiguous green region, where both classifiers simultaneously predict their respective classes with high confidence, demonstrates a fundamental limitation of the OvR approach and motivates the need for a unified multi-class classification framework that ensures unambiguous class assignments.
Figure 8.2: Illustration of the circular contradiction problem in One-vs-One classification with three classes. The plot shows three decision boundaries: R1R_{1} (solid blue line) separating C1C_{1} from C2C_{2} , R2R_{2} (dashed green line) separating C2C_{2} from C3C_{3} , and R3R_{3} (dash-dotted purple line) separating C1C_{1} from C3C_{3} . The pink shaded region highlights where the binary decisions create a circular contradiction: points in this region are classified as C2C_{2} over C1C_{1} by R1R_{1} , as C3C_{3} over C2C_{2} by R2R_{2} , but then as C1C_{1} over C3C_{3} by R3R_{3} . This cyclic inconsistency ( C1→C2→C3→C1C_{1}\rightarrow C_{2}\rightarrow C_{3}\rightarrow C_{1} ) demonstrates why decomposing multi-class problems into independent binary decisions is problematic. The perpendicular arrows on each boundary indicate the direction of classification decisions, showing how the circular contradiction arises from the interaction of these independent binary classifications.
Figure 8.3: Visualization of how softmax-based multi-class classification ensures consistent decisions without ambiguity. The left panel shows the probability surfaces P(Ck∣x)P(C_{k}|\mathbf{x}) for each class kk , demonstrating how the softmax function transforms linear predictions wkTx\mathbf{w}_{k}^{T}\mathbf{x} into proper probability distributions that sum to 1 at every point in feature space. The right panel projects these surfaces onto the feature plane, where colored regions indicate the dominant class (highest probability), solid lines mark decision boundaries where the highest probability changes from one class to another, and dashed lines show where any two classes have equal probabilities. Unlike the One-vs-One approach in Figure 2, there can be no ambiguity because every point has a well-defined probability distribution over all classes through the shared softmax denominator, with decision boundaries emerging naturally from comparing these interdependent probabilities rather than from independent binary classifications.
Figure 8.4: Extension of the softmax visualization to four classes, demonstrating how the framework naturally scales to higher dimensions. The left panel shows four probability surfaces P(Ck∣x)P(C_{k}|\mathbf{x}) that sum to 1 at every point, while the right panel reveals the resulting decision structure with (42)=6\binom{4}{2}=6 pairwise boundaries. Solid lines mark the actual decision boundaries where the highest probability changes between classes, while dashed lines indicate where pairs of classes have equal probabilities. The colored regions show where each class achieves the highest probability. Despite the increased complexity with four classes, the softmax function still ensures a consistent probability distribution everywhere in feature space, eliminating any possibility of contradictory classifications. This window-pane pattern of boundaries emerges naturally from comparing probabilities in our unified framework, rather than from independent binary decisions.
Figure 8.5: Comparison of (Pearson) correlation and mutual information for detecting statistical dependence. Left panel: Scatter plot showing data points arranged in a circular pattern. Despite the clear geometric structure, the correlation is nearly zero ( ρ≈0\rho\approx 0 ) because the relationship is not linear - as XX increases, YY sometimes increases and sometimes decreases, averaging out in the correlation calculation. Right panel: The product of marginal distributions p(x)p(y)p(x)p(y) shows what we would expect if XX and YY were independent. The significant difference between this independent model and the actual circular pattern in the data (left panel) leads to high mutual information ( I(X;Y)≈1I(X;Y)\approx 1 bit), correctly identifying the strong statistical dependence between the variables. This example illustrates how mutual information can detect nonlinear relationships that correlation misses.
Figure 9.1: Comparison between standard and Bayesian approaches to logistic regression. Left panel: The maximum likelihood solution produces a single confident decision boundary that can lead to incorrect classifications in regions with sparse or no training data. The yellow star indicates an object that is incorrectly classified with high confidence despite being in a region without training examples. Right panel: The Bayesian approach acknowledges uncertainty in the decision boundary, particularly in regions far from training data. The purple contour lines show the predictive uncertainty (standard deviation), with higher values in areas where the model should be less confident. This illustrates how Bayesian logistic regression naturally expresses uncertainty in regions of feature space that are poorly sampled by training data, avoiding overconfident but incorrect classifications for novel objects. This property is especially valuable in astronomical classification tasks where we frequently encounter objects unlike anything in our training set.
Figure 9.2: Visualization of the Laplace approximation principle. The true posterior distribution (blue solid line) is non-Gaussian with asymmetric tails, which is typical for posteriors in logistic regression with limited data. The Laplace approximation (orange dashed line) produces a Gaussian that matches the true distribution at its mode (black dot) and captures the local curvature around this peak. The approximation is determined by two key properties: (1) the location of the mode, which becomes the mean of the approximating Gaussian, and (2) the curvature at the mode (second derivative of the log posterior), which determines the variance of the Gaussian. This approach provides an excellent approximation near the mode where most probability mass is concentrated, especially as the sample size increases and the posterior becomes more concentrated.
Figure 9.3: Visualization of the Bernstein-von Mises theorem, demonstrating why Laplace approximation becomes increasingly accurate with larger sample sizes. Left panel: With only 5 observations, the posterior distribution (blue) is noticeably asymmetric, and the Gaussian approximation (orange dashed line) provides a poor fit, particularly in the tails. Right panel: With 200 observations, the posterior becomes nearly indistinguishable from a Gaussian distribution, making the Laplace approximation highly accurate. The vertical green dashed line indicates the true parameter value. This convergence to normality occurs regardless of the prior distribution, and extends to multivariate posteriors in models like logistic regression.
Figure 9.4: Illustration of how boundary points affect uncertainty in Bayesian logistic regression. Left panel: With only the original data points, there is significant uncertainty in the decision boundary, shown by the wide purple uncertainty region and the spread of possible boundaries. Right panel: After adding data points specifically near the decision boundary (triangular markers), the uncertainty region narrows significantly, and possible boundaries become more tightly clustered. This demonstrates why points near the decision boundary provide the most information about model parameters. The information value of a data point is proportional to σ(wMAPTxn)(1−σ(wTxn))\sigma(\mathbf{w}_{\text{MAP}}^{T}\mathbf{x}_{n})(1-\sigma(\mathbf{w}^{T}\mathbf{x}_{n})) , which is maximized when σ(wMAPTxn)=0.5\sigma(\mathbf{w}_{\text{MAP}}^{T}\mathbf{x}_{n})=0.5 (exactly at the boundary) and approaches zero for points far from the boundary.
Figure 9.5: Bayesian logistic regression illustrated through multiple decision boundaries sampled from the posterior distribution. Each purple line represents a possible decision boundary according to the posterior, while the bold black line shows the MAP estimate. The yellow star indicates a test point in an uncertain region, with the inset histogram displaying the distribution of predictions from all sampled models. Unlike standard logistic regression which produces a single confident boundary, the Bayesian approach naturally quantifies uncertainty through this distribution of predictions. Note how predictions span both classes (blue for Stars, red for Galaxies), reflecting genuine ambiguity in regions far from training data.
Figure 9.6: Visualization of the dimensionality reduction from weight space to logit space in Bayesian logistic regression. Panel (a) shows the 2D weight space where each point represents a possible weight vector w\mathbf{w} . The blue arrow indicates the direction of the test point x∗\mathbf{x}^{*} , while the green dashed arrow shows the direction perpendicular to x∗\mathbf{x}^{*} . The dashed lines represent iso- aa contours, where all weight vectors along a given contour yield the same value of a=wTx∗a=\mathbf{w}^{T}\mathbf{x}^{*} . Panel (b) presents the same relationship as a 3D surface, where the height corresponds to the value of aa , with red dashed lines highlighting constant- aa contours. This geometric perspective explains why we can simplify the integration from the high-dimensional weight space to a one-dimensional integral over aa : many different weight vectors map to the same scalar aa value, and for logistic regression, the prediction depends only on this scalar. This insight enables the mathematical transformation ∫σ(wTx∗)p(w∣D)dw=∫σ(a)p(a)da\int\sigma(\mathbf{w}^{T}\mathbf{x}^{*})p(\mathbf{w}|\mathcal{D})d\mathbf{w}=\int\sigma(a)p(a)da , dramatically simplifying the computation of the predictive distribution.
Figure 9.7: The probit function as the cumulative distribution function (CDF) of the standard normal distribution. Panel (a) shows the unit normal probability density function (PDF) N(x∣0,1)\mathcal{N}(x|0,1) . The shaded orange area represents the cumulative probability up to x=1x=1 , which equals Φ(1)=0.841\Phi(1)=0.841 . Panel (b) displays the probit function Φ(x)\Phi(x) , which is the CDF of the standard normal distribution. For any value xx , the probit function Φ(x)\Phi(x) gives the probability mass contained in the normal distribution up to that point.
Figure 9.8: Comparison of the sigmoid function σ(x)\sigma(x) with different scaled probit functions Φ(λx)\Phi(\lambda x) . The optimal scaling factor λ=π/8≈0.63\lambda=\sqrt{\pi/8}\approx 0.63 (red dashed line) provides the best overall approximation to the sigmoid function (blue solid line). This scaling factor is derived by matching the derivatives of both functions at x=0x=0 , where they cross the probability value of 0.5. The inset shows the absolute error for each approximation across the input range, demonstrating that the optimal scaling minimizes the maximum error. Alternative scaling factors that are either too small ( λ=0.4\lambda=0.4 , green) or too large ( λ=0.9\lambda=0.9 , orange) produce poorer approximations. This optimal approximation enables us to transform the intractable sigmoid-Gaussian integral in the predictive distribution into an analytically solvable probit-Gaussian integral while maintaining high accuracy, particularly near the decision boundary.
Figure 9.9: Comparison of standard and Bayesian logistic regression predictions. Panel (a) shows the standard logistic regression probability contours with a sharp decision boundary and increasing confidence away from the boundary regardless of data density. Panel (b) shows Bayesian logistic regression probability contours, which feature a more gradual transition between classes in regions far from training data. Panel (c) visualizes the difference between Bayesian and standard predictions, with blue regions indicating where Bayesian predicts lower probabilities (more likely to be stars) and red regions indicating higher probabilities (more likely to be galaxies) compared to standard logistic regression. The white regions in the center show where predictions are pulled toward 0.5 (maximum uncertainty). This occurs in the sparse regions between and away from the training clusters, illustrating how the uncertainty factor κ(σa2)\kappa(\sigma_{a}^{2}) naturally provides a measure of confidence that depends on the distribution of training data. The decision boundary in Bayesian logistic regression (red dashed line) is essentially identical to the standard approach (blue solid line) where data is plentiful, but the confidence in predictions away from data is appropriately reduced.
Figure 10.1: Vector projection in 2D space, illustrating a key mathematical concept for PCA. The figure shows an original vector x\mathbf{x} (blue), a unit vector b\mathbf{b} (red), and the projection x^\widehat{\mathbf{x}} of x\mathbf{x} onto b\mathbf{b} (green). The equation x^=(xTb)b\widehat{\mathbf{x}}=(\mathbf{x}^{T}\mathbf{b})\mathbf{b} shows how we compute this projection: first calculating the scalar projection xTb\mathbf{x}^{T}\mathbf{b} (which measures how much of x\mathbf{x} extends in direction b\mathbf{b} ), then multiplying by b\mathbf{b} to get a vector with the appropriate magnitude in the direction of b\mathbf{b} . The dashed line represents the perpendicular distance from x\mathbf{x} to its projection, which is minimized by this projection. This minimization property is crucial to PCA, where we seek basis vectors that minimize the sum of squared projection errors across all data points.
Figure 10.2: Visualization of Principal Component Analysis (PCA) for dimension reduction. Left panel: Original 3D data with the principal component plane. Blue points show the original data distributed around a 2D plane in 3D space. Orange arrows represent the first two principal components (PC1 and PC2), which define the optimal projection plane. Dashed lines connect original points to their projections (hollow orange circles) on the PC plane. Right panel: The reduced 2D representation of the data after projection onto the principal component plane. Points are colored by their original Z-value in 3D space, showing how the dimensionality reduction preserves the structure of the data. This illustrates how PCA identifies the lower-dimensional manifold that best preserves the information content of high-dimensional data.
Figure 10.3: Geometric intuition of Lagrange multipliers. The contour lines represent the function f(x,y)=x2+y2f(x,y)=x^{2}+y^{2} , with the constraint g(x,y)=x+y−1=0g(x,y)=x+y-1=0 shown as the colored line. Colors along the constraint indicate function values, with warmer colors representing higher values. At the constrained extremum (red dot), the gradient ∇f\nabla f (blue arrow) is perpendicular to the tangent of the constraint (green arrow). This illustrates a key insight of Lagrange multipliers: at a constrained extremum, moving along the constraint path doesn’t change the function value to first order. The figure also shows the unconstrained minimum (purple dot) at the origin, which is inaccessible due to the constraint.
Figure 10.4: Perpendicularity of the constraint gradient. This figure illustrates why the gradient of the constraint function ∇g\nabla g (blue arrows) is always perpendicular to the constraint path g(x,y)=0g(x,y)=0 (red line). The contour lines represent different values of the constraint function g(x,y)=x+y−1g(x,y)=x+y-1 , with the zero-valued contour highlighted in red. Green arrows show tangent vectors to the constraint path at various points. Because the gradient of any function points in the direction of maximum increase, and the constraint function must remain constant (equal to zero) along the constraint path, the gradient ∇g\nabla g must be perpendicular to any movement along the constraint. This property explains why, at a constrained extremum, the condition ∇f=λ∇g\nabla f=\lambda\nabla g ensures that ∇f\nabla f is perpendicular to the constraint path, which is necessary for a constrained extremum.
Figure 10.5: Geometric visualization of how eigenvectors maintain their direction under covariance transformation. The figure shows vectors on a unit circle (black dashed line) and their transformation by a covariance matrix, resulting in an ellipse (blue solid line). Eigenvectors (red and orange) preserve their directions exactly when transformed (shown in lighter shades), being scaled by their respective eigenvalues but not rotated. In contrast, non-eigenvectors (blue) change both direction and length when transformed (light blue). This property makes eigenvectors uniquely valuable for PCA: they represent invariant directions in the data space where the covariance transformation produces pure scaling without rotation, corresponding precisely to the principal axes of variation in the data.
Figure 10.6: Visualization of how PCA identifies principal components using the covariance matrix. The blue points represent the original data, and the light blue ellipse shows the covariance structure (at 2 σ\sigma ). The red and orange arrows represent the eigenvectors of the covariance matrix, with their lengths scaled proportionally to their corresponding eigenvalues ( λ1\lambda_{1} and λ2\lambda_{2} ). The first principal component (red arrow) points in the direction of maximum variance in the data, while the second principal component (orange arrow) is orthogonal to the first and captures the remaining variance. The percentages shown in the figure indicate the amount of total variance explained by each principal component. This illustrates how PCA decomposes the covariance structure of data to find an optimal lower-dimensional representation that preserves maximum variance while reducing dimensionality.
Figure 10.7: Visualization of residual data after extracting the first principal component. Blue points represent the original centered data, while red points show the residual data X^\widehat{\mathbf{X}} after projecting out the component along the first principal component b1\mathbf{b}_{1} (orange arrow). Dashed lines indicate the projection process. The semi-transparent orange plane represents the hyperplane orthogonal to b1\mathbf{b}_{1} , where all residual data points lie. This illustrates a key concept from the induction proof: after removing the first principal component, the residual data has zero variance in the b1\mathbf{b}_{1} direction. The residual covariance matrix S^\widehat{\mathbf{S}} has the first eigenvector b1\mathbf{b}_{1} with eigenvalue zero, while subsequent eigenvectors maintain their original eigenvalues. This demonstrates why the largest eigenvalue of S^\widehat{\mathbf{S}} corresponds to the second largest eigenvalue of the original covariance matrix, proving that successive principal components are eigenvectors with decreasing eigenvalues.
Figure 10.8: Scree plot for selecting the optimal number of principal components. The blue bars and orange line show the individual explained variance percentage for each principal component, demonstrating how the contribution of successive components diminishes. The green line displays the cumulative explained variance, reaching thresholds of 90%, 95%, and 99% (gray horizontal lines) at different component numbers. This visualization guides dimension selection by revealing the trade-off between dimensionality reduction and information preservation. The steep initial decline followed by flattening curve is characteristic of datasets with a few dominant modes of variation, typical in astronomical applications where high-dimensional data often arises from a smaller number of underlying physical processes.
Figure 10.9: Illustration of PCA’s limitation with multimodal data. Left panel: Bimodal data with two distinct clusters (blue and orange), each with its own internal structure and major axis (shown as dashed arrows). Notice how the first principal component (PC1, green arrow) does not align with the major direction of either cluster. Instead, PC1 points between the clusters to maximize global variance across the entire dataset. This misalignment occurs because PCA treats all data points as samples from a single distribution rather than recognizing the bimodal structure. Right panel: Data projected onto the principal components, showing how PC1 primarily captures the between-cluster variation rather than the within-cluster structure. This demonstrates why PCA can be inappropriate for astronomical datasets with natural clustering, such as distinct galaxy types or separate stellar populations. For such multimodal astronomical data, techniques that first identify clusters, as we will discuss in the next chapter, before performing dimensionality reduction often yield more physically meaningful representations.
Figure 11.1: Comparison of classification (supervised learning) and clustering (unsupervised learning) approaches in machine learning. Left panel: Classification with known labels, where data points are assigned to predefined classes (blue circles for Class A, orange squares for Class B, and green triangles for Class C). Decision boundaries (dashed lines) separate the feature space into regions for different classes based on labeled training data. The horizontal line separates Class C from Classes A and B, while the vertical line separates Class A from Class B. Right panel: Clustering without prior labels, where algorithms discover natural groupings in unlabeled data (gray circles). Discovered cluster centers (colored X marks) and their corresponding regions (colored ellipses) emerge from the data’s inherent structure rather than predefined categories.
Figure 11.2: Visualization of the dual optimization challenge in K-means clustering. The algorithm simultaneously optimizes two sets of parameters: (1) the cluster centers μk\boldsymbol{\mu}_{k} (marked with X symbols) and (2) the binary assignment variables rik∈{0,1}r_{ik}\in\{0,1\} that assign each data point to exactly one cluster. The objective function J=∑i=1N∑k=1Krik∥xi−μk∥2J=\sum_{i=1}^{N}\sum_{k=1}^{K}r_{ik}\|\mathbf{x}_{i}-\boldsymbol{\mu}_{k}\|^{2} represents the sum of squared distances between data points and their assigned centers, which the algorithm aims to minimize. Dashed lines connect selected data points to their assigned cluster centers, with rik=1r_{ik}=1 indicating that point ii belongs to cluster kk . This presents a circular dependency in the optimization process: the optimal centers depend on the point assignments, while the optimal assignments depend on the center positions. This interdependence necessitates the iterative Expectation-Maximization approach used in K-means, where we alternate between updating assignments and center positions until convergence.
Figure 11.3: Visualization of the Expectation-Maximization (EM) algorithm for K-means clustering across three iterations. Top row: Initial centers (left) are chosen arbitrarily, then the first E-step (right) assigns each point to its nearest center (shown by colors and dashed lines). Middle row: The first M-step (left) moves each center to the mean of its assigned points, with text showing the exact coordinate transition. The second E-step (right) then reassigns points based on these updated centers. Bottom row: The second M-step (left) further refines the centers, followed by the third E-step (right) which shows the final assignments. This demonstrates how K-means resolves the circular dependency between cluster centers and point assignments through iterative optimization.
Figure 11.4: Demonstration of K-means sensitivity to initialization. Left: Good initialization with centers placed near each true cluster leads to correct identification of all three natural clusters. Initial centers (circles) converge to optimal final positions (X markers). Right: Poor initialization with two centers in the top cluster and one center between the bottom clusters results in a suboptimal solution where the large top cluster is incorrectly split while the two smaller bottom clusters are incorrectly merged. This illustrates how K-means can converge to local minima when initial centers are poorly positioned, highlighting why strategies like multiple random initializations or K-means++ are often necessary in practice.
Figure 11.5: Illustration of the K-means++ initialization algorithm. Step 1 (Top Left): The first center is chosen uniformly at random from the data points. Step 2 (Top Right): The second center is selected with probability proportional to the squared distance from the first center, visualized by circle size. This gives distant points higher probability of selection, encouraging diversity in center placement. Step 3 (Bottom Left): The third center is selected with probability proportional to the squared distance from the nearest existing center, again using weighted sampling. Step 4 (Bottom Right): The final initialization shows the three centers with their corresponding assignment regions. K-means++ addresses the initialization sensitivity in standard K-means by maximizing the minimum distance between centers, leading to faster convergence and better clustering results compared to random initialization, particularly when distinct cluster structures exist in the data.
Figure 11.6: Visualization of the elbow method for determining the optimal number of clusters (K) in K-means. Top Left: The original dataset with four natural clusters. Top Right: K-means clustering results with K=4, showing the identified clusters and their centers (red X markers). Bottom Left: The elbow plot showing inertia (within-cluster sum of squares) versus K. The characteristic “elbow” at K=4 indicates the point where adding more clusters yields diminishing returns. This is analogous to the scree plot in PCA for determining the number of components to retain. Bottom Right: Additional analysis showing the rate of change of inertia (approximate derivative) and silhouette scores across different K values, both confirming K=4 as optimal. The silhouette score peaks at K=4, indicating that points are well-matched to their own clusters and well-separated from other clusters at this value.
Figure 11.7: Demonstration of K-means clustering limitation with non-convex data structures. The data consists of two half-moon shaped clusters (distinguished by different marker shapes: circles and triangles), but K-means with K=2 creates an inappropriate linear decision boundary (dashed line) that fails to properly separate these naturally curved structures. The colors represent K-means cluster assignments, which clearly cut across the true half-moon patterns rather than respecting their intrinsic shapes. This limitation stems directly from the spherical and convex cluster assumption inherent in K-means, where decision boundaries are always hyperplanes perpendicular to the line connecting cluster centers, similar to logistic regression.
Figure 11.8: Illustration of Gaussian Mixture Models (GMM) with three clusters demonstrating soft assignment. Unlike K-means’ hard boundaries, GMM assigns each point a probability of belonging to each cluster. Points are colored according to their most likely cluster assignment, with cluster centers marked by X symbols. The elliptical contours around each center show the 1, 2, and 3 standard deviation boundaries of each Gaussian component, highlighting how GMMs can model clusters with different shapes, sizes, and orientations through their covariance matrices. Red-outlined points indicate ambiguous cases where no single cluster has a dominant ( >80%>80\% ) probability of membership. This probabilistic framework is a key advantage of GMMs over K-means, allowing them to capture uncertainty in cluster assignments and better represent the underlying data distribution, particularly in regions where clusters overlap.
Figure 11.9: Illustration of mixture weights ( πk\pi_{k} ) in Gaussian Mixture Models. A 1D dataset generated from three Gaussian components with different weights ( π1=0.50\pi_{1}=0.50 , π2=0.31\pi_{2}=0.31 , π3=0.19\pi_{3}=0.19 ) is shown as a histogram. Each colored curve represents a component density scaled by its corresponding mixture weight. The mixture weights represent the prior probability that a data point belongs to each component before observing its features, and they determine the relative contribution of each Gaussian to the overall mixture. These weights must satisfy πk≥0\pi_{k}\geq 0 for all kk and ∑k=1Kπk=1\sum_{k=1}^{K}\pi_{k}=1 , forming a proper probability distribution over the components.
Figure 11.10: Comparison of covariance constraints in Gaussian Mixture Models using the same dataset with correlated features. Left: Full covariance matrices ( Σk\boldsymbol{\Sigma}_{k} ) with K⋅D(D+1)/2K\cdot D(D+1)/2 parameters can capture correlations between features, allowing ellipses of any orientation to accurately model the true cluster shapes. Middle: Diagonal covariance matrices ( Σk=diag(σk12,σk22,…,σkD2)\boldsymbol{\Sigma}_{k}=\text{diag}(\sigma_{k1}^{2},\sigma_{k2}^{2},\ldots,\sigma_{kD}^{2}) ) with K⋅DK\cdot D parameters constrain ellipses to align with coordinate axes, preventing the model from capturing feature correlations. This represents scenarios where different features vary independently within each cluster. Right: Spherical covariance matrices ( Σk=σk2I\boldsymbol{\Sigma}_{k}=\sigma_{k}^{2}\mathbf{I} ) with just KK parameters constrain clusters to be circular, similar to K-means but maintaining the probabilistic framework. These constraints illustrate the trade-off between model flexibility and computational efficiency.
Figure 11.11: Expectation-Maximization (EM) algorithm for Gaussian Mixture Models. Top row: Initial parameters with spherical covariances (left) and the first E-step showing responsibilities γik\gamma_{ik} as color intensities (right). Middle row: First M-step updating means, covariance matrices, and weights (left), followed by the second E-step with refined responsibilities (right). Bottom row: Second M-step further refining parameters (left) and final E-step showing a point with membership in multiple components highlighted in red (right). Unlike K-means’ hard assignments, GMM uses soft probabilistic assignment through the responsibility values γik\gamma_{ik} , capturing clusters with different shapes and orientations through component-specific covariance matrices.
Figure 11.12: Comparison of Gaussian Mixture Model results with different initializations. Initial component centers are shown as hollow circles, final centers as X markers, with arrows indicating movement during EM optimization. Left: Poor initialization leading to a suboptimal fit where the algorithm fails to properly identify the true underlying clusters. Right: Better initialization resulting in a fit that accurately captures the true cluster structure. The covariance ellipses show the shape and orientation of each Gaussian component. This comparison demonstrates how the EM algorithm for GMMs can converge to different local optima depending on initialization, highlighting the importance of proper starting conditions in mixture modeling.
Figure 11.13: Using Gaussian Mixture Models for outlier detection. Left: A GMM with two components fitted to data (gray points) with outliers (red X markers) manually added. Ellipses show the 1, 2, and 3 standard deviation contours for each Gaussian component. Right: Log-likelihood landscape visualized as a heatmap, where points with low likelihood values (darker colors) are potential outliers. The dashed red contour represents an outlier threshold at the 5th percentile of likelihood values. This demonstrates how the probability density p(x)p(\mathbf{x}) provided by GMMs can naturally identify data points that lie in low-density regions and are unlikely to have been generated by the estimated model, making GMMs valuable for anomaly detection in astronomical data where outliers often represent scientifically interesting objects.
Figure 11.14: Model selection for Gaussian Mixture Models using information criteria. The plot shows how AIC (blue) and BIC (orange) values change as a function of the number of components (K). Both criteria balance the model fit (log-likelihood) against model complexity (number of parameters). The BIC applies a stronger penalty term Pln⁡(N)P\ln(N) compared to AIC’s penalty of 2P2P , where PP is the number of parameters and NN is the sample size. Consequently, BIC typically favors simpler models (lower K) than AIC. AIC tends to identify more clusters because it focuses on optimizing prediction accuracy rather than finding the “true” model, making it more willing to accept additional complexity when it offers even modest improvements in fit. The minimum of each curve indicates the optimal number of components according to that criterion. After reaching their minima, both criteria begin to increase as the penalty for additional complexity outweighs the diminishing improvements in fit. This demonstrates how information criteria implement Occam’s razor by mathematically formalizing the trade-off between model fit and complexity.
Figure 12.1: Illustration of Buffon’s needle experiment, an early example of Monte Carlo methods from the 18th century. The experiment involves randomly dropping needles of length LL (shown in vertical orientation) on a surface with parallel lines spaced at distance DD . Needles are colored according to whether they cross a line (red) or not (blue). The probability of a needle crossing a line is directly related to π\pi through the formula π≈2LN/(Dn)\pi\approx 2LN/(Dn) , where NN is the total number of needles dropped and nn is the number of needles that cross a line. This physical randomization process provides a method to empirically estimate π\pi without requiring complex mathematical calculations, illustrating the core principle behind modern Monte Carlo integration techniques.
Figure 12.2: Illustration of grid-based sampling for a bimodal Gaussian mixture. Top panel: Discretization of the continuous PDF p(x)=0.7⋅N(x∣−2,0.52)+0.3⋅N(x∣2,12)p(x)=0.7\cdot\mathcal{N}(x|-2,0.5^{2})+0.3\cdot\mathcal{N}(x|2,1^{2}) into equally spaced bins of width Δx\Delta x , with each bin shown in a different color. The black curve represents the true continuous PDF. Bottom panel: The stacked bar representation of bin probabilities, where the width of each colored segment corresponds to the probability mass of the matching bin in the top panel. The thresholds sis_{i} partition the unit interval according to these probabilities. When sampling, a uniform random number u∼U(0,1)u\sim\mathcal{U}(0,1) maps to the xx value from the bin it falls within, as demonstrated by the three colored arrows. This grid-based approach allows direct sampling from complex distributions through discretization and inverse transform sampling.
Figure 12.3: Visualization of the inverse CDF method for sampling from a bimodal Gaussian mixture. Top panel: The probability density function (PDF) p(x)=0.7⋅N(x∣−2,0.52)+0.3⋅N(x∣2,12)p(x)=0.7\cdot\mathcal{N}(x|-2,0.5^{2})+0.3\cdot\mathcal{N}(x|2,1^{2}) , with reference points highlighted in different colors. The navy shaded area represents the total probability mass. Bottom panel: The cumulative distribution function (CDF) F(x)=∫−∞xp(t)dtF(x)=\int_{-\infty}^{x}p(t)dt and illustration of the inverse CDF sampling process. The sampling procedure is demonstrated with three examples: starting with a uniform random number u∼U(0,1)u\sim\mathcal{U}(0,1) on the y-axis, tracing horizontally to the CDF curve, and projecting vertically to the x-axis to obtain the corresponding sample value. This method transforms uniform random numbers into samples that follow the target distribution by leveraging the CDF’s property of rising more steeply in high-probability regions and more gradually in low-probability regions.
Figure 12.4: Illustration of the rejection sampling method for a bimodal Gaussian mixture distribution. Top panel: The target distribution (black curve with navy shading), the proposal distribution q(x)q(x) (green dashed line), and the scaled proposal M⋅q(x)M\cdot q(x) (solid green line). The red shaded area represents the rejection region—points that fall between the target and scaled proposal distributions. Bottom panel: Demonstration of the accept/reject mechanism. For each candidate point sampled from the proposal distribution, its vertical position is determined by a uniform random number u∼U(0,1)u\sim\mathcal{U}(0,1) scaled by M⋅q(x)M\cdot q(x) . Points falling below the target density curve (teal circles) are accepted, while those in the rejection region (orange X marks) are rejected. The scatter plot at the bottom shows the distribution of the resulting samples, illustrating how rejection sampling effectively “sculpts” the proposal distribution to match the shape of the target distribution.
Figure 12.5: Importance sampling example. Top panel: The target distribution (truncated normal for x≥5x\geq 5 ) and two proposal distributions: a standard normal distribution (Proposal 1) and a shifted exponential distribution (Proposal 2). The standard normal proposal places very little probability mass in the region of interest, resulting in a very low effective sample size (ESS ≈3\approx 3 for N=107N=10^{7} samples). In contrast, the shifted exponential proposal concentrates its probability mass in the region of interest, leading to a much higher effective sample size (ESS ≈0.8N\approx 0.8N ). Bottom panel: Convergence of importance sampling estimates to the true mean as a function of sample size. The shifted exponential proposal (Proposal 2) converges quickly with relatively few samples, while the standard normal proposal (Proposal 1) requires billions of samples to achieve comparable accuracy. The ESS directly impacts the variance of our estimates, with the standard error scaling as O(1/ESS)\mathcal{O}(1/\sqrt{\text{ESS}}) .
Figure 13.1: Illustration of Markov Chain Monte Carlo (MCMC) sampling on a multimodal posterior distribution. The background contours show the probability density of a complex posterior with three distinct modes. The orange line represents the path of the MCMC chain as it explores the parameter space, with small directional arrows indicating the progression of the chain. The chain continuously navigates through the parameter space, exploring both high and low probability regions, though it naturally spends more time in areas of higher probability density. This sequential exploration strategy transforms an intractable sampling problem into a tractable guided random walk, allowing MCMC to efficiently sample from complex distributions that would be impossible to characterize with simpler methods.
Figure 13.2: Visualization of the transition dynamics in Markov Chain Monte Carlo. The blue contours represent the posterior distribution p(x∣D)p(\mathbf{x}|\mathcal{D}) that we aim to sample from. The current state (dark blue point) is the center of a transition matrix (green dashed contours) that proposes possible next states. Proposals toward regions of higher posterior density (light blue points) have high acceptance probabilities, indicated by thicker arrows, while proposals toward regions of lower posterior density (orange points) have low acceptance probabilities, shown with thinner arrows. This mechanism ensures that the Markov chain spends more time in high-probability regions while occasionally exploring less probable areas, analogous to our city traveler spending more time in densely populated neighborhoods. Over many iterations, this exploration strategy generates samples that represent the target posterior distribution.
Figure 13.3: Illustration of convergence to the stationary distribution in MCMC. The blue solid curve represents the target stationary distribution, while the orange dashed line shows the initial distribution. The blue dash-dot line depicts an intermediate state during convergence. Before reaching stationarity, there exists a net flow of probability mass toward the target distribution—similar to population migration in our city analogy, with more people moving into densely populated neighborhoods than leaving them. This net flow is driven by asymmetric transition probabilities: transitions from low-probability to high-probability regions (point B to A, thick red arrow) occur with higher probability than transitions from high-probability to low-probability regions (point A to B, thin orange arrow). Once the chain reaches the stationary distribution, these flows perfectly balance each other—the detailed balance principle ensures that for every pair of states, the probability flow in each direction is exactly equal, maintaining the distribution in perfect equilibrium with no further net change. The Markov chain thus first evolves toward the target distribution and then remains there indefinitely.
Figure 13.4: Visualization of detailed balance using a discrete state space representation. Each histogram bin represents a state in the Markov chain, with heights proportional to their stationary probabilities p(i)p(i) . The arrows between states 3 and 6 show transitions in both directions, with their respective transition probabilities T(i→j)T(i\rightarrow j) and T(j→i)T(j\rightarrow i) . Detailed balance ensures that the flow from state ii to state jj equals the flow from jj to ii : p(i)⋅T(i→j)=p(j)⋅T(j→i)p(i)\cdot T(i\rightarrow j)=p(j)\cdot T(j\rightarrow i) . This equality is demonstrated numerically in the figure. When detailed balance holds for all pairs of states, the distribution remains stationary—there is no net change in the probability mass. Moreover, if the system is out of equilibrium, detailed balance actively drives it toward the equilibrium state.
Figure 13.5: Visualization of irreducibility in Markov chains. The left panel shows an irreducible chain where any state can be reached from any other state in a finite number of steps (highlighted path shows one possible route from State 1 to State 4). This property ensures that the chain can explore the entire state space, regardless of the starting point. The right panel illustrates a reducible chain where the state space is disconnected into two separate components: States 1-3 (blue) and States 4-5 (orange). States in one component cannot reach states in the other component, violating irreducibility. In astronomical applications, this would correspond to a posterior distribution with completely separated modes with zero probability between them—for example, two solutions that are equally valid but with no connecting path through parameter space. Such situations violate ergodicity and prevent conventional MCMC methods from correctly sampling the full posterior distribution.
Figure 13.6: Comparison of periodic and aperiodic Markov chains. The left panel shows a periodic chain with period 3, where states must follow a fixed cycle (1→2→3→1…). The sequence below demonstrates how the visitation pattern repeats exactly every three steps, preventing convergence to a stationary distribution. The right panel shows how adding a single self-loop to state 2 breaks this periodicity. Aperiodicity ensures that the chain doesn’t get trapped in oscillating patterns, establishing a stable visitation pattern that matches the target distribution.
Figure 13.7: Illustration of positive recurrence in Markov chains. The left panel shows a positively recurrent chain in a bounded parameter space where the chain periodically returns to previously visited states (red circles marking returns to the starting state). This property ensures that the expected time to return to any state is finite, allowing the chain to establish a proper stationary distribution. The right panel demonstrates a non-positively recurrent chain in an unbounded parameter space, where the chain can drift away from its starting point. It’s important to note that an unbounded space doesn’t necessarily imply non-positive recurrence—this is a property of the transition matrix itself, not just the parameter space. Without positive recurrence, a Markov chain cannot establish a stable equilibrium distribution, making it impossible to reliably sample from the target posterior distribution.
Figure 13.8: Step-by-step illustration of the Metropolis algorithm sampling from a multimodal distribution. Step 1 (Initialize): We choose a starting point x(0)\mathbf{x}^{(0)} in parameter space. Step 2 (Propose): A candidate state x′\mathbf{x}^{\prime} is generated from a symmetric proposal distribution (green dashed contours) centered at the current state. Step 3 (Compute): We calculate the acceptance probability as the ratio of densities, A(x(t−1),x′)=min⁡(1,p(x′)p(x(t−1)))A(\mathbf{x}^{(t-1)},\mathbf{x}^{\prime})=\min\left(1,\frac{p(\mathbf{x}^{\prime})}{p(\mathbf{x}^{(t-1)})}\right) . Step 4 (Accept/Reject): We draw a uniform random number u∼U(0,1)u\sim\mathcal{U}(0,1) and accept the proposal if u≤Au\leq A ; otherwise, we reject it and stay at the current position. We keep repeating this process for many iterations, which allows the algorithm to explore the parameter space by preferentially moving toward higher-probability regions while occasionally accepting moves to lower-probability regions, ultimately sampling from the target distribution.
Figure 13.9: Visualization of the burn-in concept in MCMC sampling. The chain starts from a low-probability region (orange point) and initially explores the parameter space inefficiently (orange trajectory). This initial phase constitutes the burn-in period, during which samples are discarded because they reflect the arbitrary starting position rather than the target distribution. After the burn-in period, the chain settles into its stationary distribution and efficiently explores the high-probability modes of the distribution (blue trajectory), providing representative samples for inference. The colorbar shows the progression of iterations with the burn-in cutoff indicated by the horizontal dashed line.
Figure 13.10: Comparison of MCMC performance with different proposal mechanisms. Left panel: A narrow proposal distribution results in a high acceptance rate, as most proposed moves stay within regions of similar posterior probability. Right panel: A wide proposal distribution yields a low acceptance rate, as proposals frequently land in low-probability regions and are rejected. The blue contours represent the target posterior distribution, while dashed contours show the proposal distributions centered at the current state (blue dot). This visualization demonstrates the trade-off in MCMC sampling between small steps (which are readily accepted but explore slowly) and large steps (which explore faster but are more frequently rejected).
Figure 13.11: Visualization of autocorrelation in MCMC sampling. Top panel: The MCMC chain trace showing a reference point at iteration t=100t=100 (purple) and corresponding points at lags k=1k=1 , k=5k=5 , and k=20k=20 (different shades of purple). Middle panel: Scatter plots showing the relationship between states xtx_{t} and xt+kx_{t+k} for different lags. The increasing scatter away from the diagonal as lag increases illustrates decreasing correlation between samples separated by larger intervals. Bottom panel: The autocorrelation function R(k)R(k) showing how correlation decays with increasing lag. This progressive decorrelation allows us to treat different points in the chain as approximately independent samples from the target distribution once sufficient lag time has passed, a consequence of the chain’s ergodic properties.
Figure 13.12: Illustration of MCMC thinning based on proposal width and autocorrelation. Left panel: Early iterations of MCMC chains with three different proposal widths (vertically shifted for clarity). The small width chain (blue) makes tiny steps with high acceptance rate but explores the parameter space very slowly. The medium width chain (red) and large width chain (green) show markers indicating which samples would be kept after appropriate thinning, with thinning factors determined by their autocorrelation properties. Right panel: Corresponding autocorrelation functions with vertical dashed lines indicating the thinning factors for the medium and large width chains, chosen where autocorrelation drops below 0.1 (horizontal dashed line). Conversely, an overly large proposal width means that most of the proposed steps are rejected during the chain, resulting in many duplicate samples where the chain remains stuck at the same position, which also reduces efficiency and effective sample size.
Figure 13.13: Illustration of the Gelman-Rubin diagnostic components and their evolution in MCMC sampling. Top panel: Four chains starting from different widely separated points gradually converge toward the target distribution. The small proposal width deliberately slows convergence to clearly demonstrate the diagnostic process. Bottom panel: Evolution of the within-chain variance ( WW ) and between-chain variance ( BB ) calculated over sliding windows. Initially, BB is much larger than WW , as chains are exploring different regions of the parameter space. As sampling progresses, BB decreases substantially while WW stabilizes, indicating chains are mixing across the same regions. Convergence is achieved when BB approaches WW (dashed reference line), corresponding to chains that are statistically indistinguishable from each other. This visualization demonstrates why the Gelman-Rubin diagnostic is effective at assessing MCMC convergence by comparing variability within chains to variability between chains, without requiring knowledge of the true target distribution.
Figure 13.14: Illustration of the Geweke convergence diagnostic for a single MCMC chain. Top panel: A single MCMC chain starting from a point far from the target distribution. The gray shaded area represents a candidate burn-in period, with segments A (first 10% of post-burn-in samples, red) and B (last 50% of post-burn-in samples, green) used for the Geweke test. If the chain has converged, these two segments should have statistically similar means. Bottom panel: Geweke Z-scores computed for different burn-in sizes. The Z-score measures the difference between the means of segments A and B, properly adjusted for autocorrelation. Values within the gray band ( ∣Z∣<2|Z|<2 ) suggest convergence, while values outside indicate that more burn-in may be needed. This diagnostic is particularly valuable when computational constraints limit the ability to run multiple chains for the Gelman-Rubin diagnostic, allowing convergence assessment using a single chain while properly accounting for autocorrelation in the samples.
Figure 13.15: Comparison of the Metropolis algorithm (left) versus the Metropolis-Hastings algorithm (right). Both panels show the same bimodal posterior distribution (blue contours) and start from the same current state (blue point). The Metropolis algorithm uses a symmetric proposal distribution (green dashed contours), where the probability of proposing a move from x\mathbf{x} to x′\mathbf{x}^{\prime} equals the probability of the reverse move. In contrast, the Metropolis-Hastings algorithm uses a highly asymmetric proposal distribution (orange dashed contours), where q(x′∣x)≠q(x∣x′)q(\mathbf{x}^{\prime}|\mathbf{x})\neq q(\mathbf{x}|\mathbf{x}^{\prime}) . The Metropolis-Hastings algorithm compensates for the proposal asymmetry through the additional correction term q(x∣x′)q(x′∣x)\frac{q(\mathbf{x}|\mathbf{x}^{\prime})}{q(\mathbf{x}^{\prime}|\mathbf{x})} in the acceptance ratio, maintaining detailed balance despite the directional bias in the proposal mechanism. This generalization allows for more efficient exploration of complex posterior geometries while preserving the theoretical guarantees of convergence to the target distribution, assuming ergodic conditions are met.
Figure 13.16: Illustration of Gibbs sampling. The main panel shows the joint posterior distribution (blue contours) with the characteristic stair-step pattern of Gibbs sampling. Purple dashed lines indicate the conditional slices where we sample one parameter given the other. Each iteration consists of two steps: first sampling x1x_{1} given x2x_{2} (horizontal movement), then sampling x2x_{2} given the new value of x1x_{1} (vertical movement). The top and right panels show the one-dimensional conditional distributions p(x1∣x2)p(x_{1}|x_{2}) and p(x2∣x1)p(x_{2}|x_{1}) , with purple solid lines marking the modes of these conditionals and green dashed lines showing the actual values sampled. This visual demonstrates how Gibbs sampling decomposes the high-dimensional sampling problem into a series of conditional sampling tasks.
Figure 13.17: Effect of parameter correlation on Gibbs sampling efficiency. Top row: Sampling paths in parameter space for bivariate Gaussian distributions with no correlation ( ρ=0\rho=0 , left) and high correlation ( ρ=0.95\rho=0.95 , right). The characteristic rectangular pattern of Gibbs sampling becomes increasingly narrow with higher correlation, showing the algorithm’s difficulty in exploring the major axis of the distribution. Middle row: MCMC chains showing stronger autocorrelation in the high-correlation case. Note the dramatically reduced effective sample size (ESS) percentage when correlation is high. Bottom row: Autocorrelation functions for both parameters, demonstrating how autocorrelation decays much more slowly with higher parameter correlation. This example illustrates why Gibbs sampling, despite its rejection-free nature, can be highly inefficient for strongly correlated parameters. While the acceptance rate remains perfect in both cases, the high-correlation scenario requires many more iterations to effectively explore the posterior distribution.
Figure 14.1: Comparison of feature-based and kernel-based approaches for regression tasks. The top panel illustrates the feature-based approach, where inputs xx are explicitly mapped to a feature space ϕ(x)\boldsymbol{\phi}(x) before applying linear regression. The black points represent training data, while the blue dashed and orange solid lines show predictions from low-degree (underfitting) and high-degree (more flexible) polynomial models, respectively. The bottom panel demonstrates the kernel-based approach using Gaussian Processes, which implicitly works in an infinite-dimensional feature space through the kernel trick. The blue line shows the GP mean prediction with the light blue region indicating the 95% confidence interval, naturally providing uncertainty quantification. The green dashed line (right y-axis) visualizes the kernel similarity function from the reference point x0=0x_{0}=0 , illustrating how Gaussian Processes measure similarity between inputs rather than explicitly defining features.
Figure 14.2: Kernel Density Estimation (KDE) illustrating the core concept behind kernel methods. The figure shows how KDE constructs a flexible, non-parametric probability density estimate by placing a kernel function at each observed data point (black tick marks along the x-axis). The overall density estimate is created by summing these kernels, allowing complex functions to emerge from simple building blocks without requiring explicit feature engineering. The bandwidth parameter hh controls the smoothness of the estimate, showcasing the bias-variance tradeoff: the red curve (h=0.05) is undersmoothed with high variance but captures fine details, the blue curve (h=0.15) provides optimal smoothing that closely follows the true density, while the green (h=0.6) and purple (h=1.5) curves progressively oversmooth the data. This same principle—using kernels to measure similarity between points—underlies Gaussian Process regression, where inputs with similar features produce similar outputs without requiring explicit feature mappings.
Figure 14.3: Illustration of kernel-based regression through weighted averaging. This figure demonstrates how predictions can be made at a new point x∗x_{*} (red star) using weighted combinations of existing observations (black dots) without requiring an explicit parametric model. The weights are determined by a kernel function (light blue curve) centered at the prediction point, which assigns higher weights to temporally proximate observations. The thickness of the blue arrows represents the relative contribution (weight) of each observation to the final prediction, with nearby points having stronger influence. Instead of engineering explicit features or constructing a complex physical model, we define similarity between points via the kernel function and let the data speak for itself. The prediction is a weighted average y(x∗)=∑i=1Nwiyi∑i=1Nwiy(x_{*})=\frac{\sum_{i=1}^{N}w_{i}y_{i}}{\sum_{i=1}^{N}w_{i}} , where weights naturally decrease with distance from the prediction point.
Figure 14.4: Comparison of RBF kernel (left) and Matérn kernels with different smoothness parameters (right). Both panels show kernel functions k(x,x0)k(x,x_{0}) centered at the reference point x0=0x_{0}=0 . The RBF kernel exhibits infinite smoothness, decreasing exponentially with squared distance. Matérn kernels demonstrate how the smoothness parameter ν\nu controls function behavior: ν=0.5\nu=0.5 (red) produces continuous but not differentiable functions, ν=1.5\nu=1.5 (green) yields once-differentiable functions, and ν=2.5\nu=2.5 (purple) produces twice-differentiable functions. As ν\nu increases, the Matérn kernel converges to the RBF kernel. This flexibility in controlling function smoothness makes these kernels valuable for modeling astronomical phenomena where different physical processes may exhibit varying degrees of smoothness.
Figure 14.5: Gaussian Process regression demonstrating the key features of the full Bayesian posterior. The figure shows observations (black points), the GP mean prediction (blue line), and 95% confidence intervals (light blue region), along with the true underlying function (dashed gray). Note how the posterior uncertainty varies with data density: it is minimal near observed points and grows in regions where data is sparse or absent. This behavior emerges naturally from the mathematics of GP regression through the posterior variance σf2=k(x∗,x∗)−k∗T(K+σn2I)−1k∗\sigma_{f}^{2}=k(x_{*},x_{*})-\mathbf{k}_{*}^{T}(\mathbf{K}+\sigma_{n}^{2}\mathbf{I})^{-1}\mathbf{k}_{*} , which decreases when the test point x∗x_{*} is similar to many training points. This automatic uncertainty quantification makes GPs particularly valuable for astronomical applications with irregular sampling and data gaps, providing principled uncertainty estimates that prevent overconfident conclusions in sparsely observed regions.
Figure 14.6: Illustration of the difference between latent function values ( f∗f_{*} ) and observed outputs ( y∗y_{*} ) in Gaussian Process regression. The plot shows training data with error bars representing observation noise (black points), the true function (dashed gray), and GP mean prediction (blue). Two types of uncertainty bands are displayed: the inner blue band shows the 95% confidence interval for the latent function f∗f_{*} (uncertainty in the underlying pattern), while the outer red band shows the 95% confidence interval for noisy observations y∗y_{*} (includes both function uncertainty and observation noise). The mathematical relationship between these uncertainties is straightforward: Var(y∗)=Var(f∗)+σn2\text{Var}(y_{*})=\text{Var}(f_{*})+\sigma_{n}^{2} , where σn2\sigma_{n}^{2} is the noise variance. When analyzing phenomena like stellar variability or quasar light curves, astronomers are often interested in the latent function f∗f_{*} to understand the underlying physical processes, while accounting for the total uncertainty including noise when making predictions about future observations y∗y_{*} .
Figure 14.7: Visualization of function spaces defined by Gaussian Processes with different RBF kernel parameters. Each panel shows five random functions sampled from a GP prior with zero mean and a specified RBF kernel. Left: Functions sampled from a GP with longer length scale ( ℓ=1.0\ell=1.0 ) exhibit smoother behavior with gradual variations. Right: Functions sampled from a GP with shorter length scale ( ℓ=0.3\ell=0.3 ) show more rapid fluctuations and can capture finer details. Both GPs use the same amplitude parameter ( σf=1.0\sigma_{f}=1.0 ) and the same random seeds, isolating the effect of the length scale. This illustrates a key property of the function-space view of Gaussian Processes: the kernel directly determines the properties of functions in the prior distribution.
Figure 14.8: Comparison of Gaussian Process prior and posterior, illustrating the connection to conditioning multivariate Gaussians. Left panel: The GP prior shows five random functions sampled from a zero-mean GP with an RBF kernel. This is analogous to sampling from the joint distribution of a high-dimensional Gaussian before conditioning. Right panel: The GP posterior after conditioning on five observations (black points) demonstrates how observing data updates our beliefs. The posterior mean (solid blue line) and 95%95\% confidence interval (light blue region) show reduced uncertainty near observations. The inset plot directly illustrates this conditioning process through a 2D Gaussian representing the joint prior distribution of function values at two arbitrary points (xa,xb)(x_{a},x_{b}) . The red dashed line shows conditioning on f(xa)=1.0f(x_{a})=1.0 , while the red curve displays the resulting conditional distribution p(f(xb)∣f(xa))p(f(x_{b})|f(x_{a})) - a direct visualization of how conditioning one point in a GP affects correlated points based on their distance in correlation space. This reinforces that GP regression applies Gaussian conditioning principles in infinite dimensions.
Figure 14.9: Visualization of how different kernel parameters affect functions sampled from Gaussian Process priors. The figure shows random functions drawn from GP priors with varying kernel parameters. Top row: Functions sampled with smaller parameter values; bottom row: Functions with larger parameter values. From left to right: (1) RBF amplitude parameter θ0\theta_{0} controls the vertical scale of functions, with larger values allowing greater deviations from the mean; (2) Inverse length scale parameter θ1\theta_{1} determines function smoothness, with smaller values producing smoother functions and larger values allowing more rapid variations; (3) Period parameter θ4\theta_{4} in the periodic kernel sets the distance between repetitions in the function; (4) Flexibility parameter θ5\theta_{5} controls how strictly the periodicity is enforced, with smaller values requiring more exact repetition and larger values allowing the periodic pattern to evolve. Each subplot shows three random function samples (in different colors) from the same GP prior to illustrate the characteristic properties induced by each parameter setting. This visualization demonstrates how kernel parameters directly encode prior beliefs about function properties in the GP framework.
Figure 14.10: Demonstration of hyperparameter optimization in Gaussian Process regression. The top panel shows the log marginal likelihood landscape over a grid of RBF kernel hyperparameters (length scale and amplitude) plotted on logarithmic scales. Brighter regions indicate higher likelihood values, with the optimal hyperparameters (red star) yielding significantly higher likelihood than suboptimal choices (red x). The bottom panel compares GP regression results using optimal versus suboptimal hyperparameters. The optimal model (blue curve) with appropriate length scale and amplitude captures the true underlying function (dashed gray) more accurately, with well-calibrated uncertainty bands (light blue region). In contrast, the suboptimal model (orange) with shorter length scale oscillates unnecessarily between data points, exhibiting overfitting behavior. This figure illustrates how maximizing the marginal likelihood automatically balances model fit and complexity without requiring explicit regularization terms.
Figure 14.11: Comparison of decision boundaries and predictive uncertainty between Bayesian Linear Classification (left) and Gaussian Process Classification (right) on a synthetic moon-shaped dataset. Colors represent the predicted probability of class membership, with red indicating high probability of class 1 and blue indicating high probability of class 0. The solid black line shows the decision boundary where the class probability equals 0.5. The Bayesian Linear Classifier is limited by its linear nature and cannot capture the moon-shaped structure of the data. In contrast, the Gaussian Process Classifier exhibits both increased uncertainty in sparse data regions and a flexible, nonlinear decision boundary that properly adapts to the data’s intrinsic geometry. This comparison illustrates how GPC extends logistic regression beyond linear boundaries while maintaining proper uncertainty quantification—the model becomes appropriately less confident in regions with sparse or no data.
Figure 14.12: Visualization of the information flow in Gaussian Process Classification. The top panels show the progression from training data (left) through latent function values (middle) to class probabilities (right) using moon-shaped data, illustrating how GPC can handle nonlinear decision boundaries. The middle panel demonstrates how uncertainty is propagated from latent space to probability space: the posterior standard deviation in latent space ( σ∗\sigma_{*} ) represents our uncertainty about the true latent function value a∗a_{*} at a test point. When transformed through the sigmoid function, this uncertainty becomes asymmetric in probability space - wider near the decision boundary (where the latent function a∗≈0a_{*}\approx 0 ) and narrower at the extremes (where ∣a∗∣|a_{*}| is large). The bottom diagram shows the core conceptual flow: the Gaussian Process transforms input variables to latent variables (assigning a latent value aia_{i} to each input xi\mathbf{x}_{i} ), then the sigmoid response function converts these latent values to class probabilities. This structure illustrates how GPC extends logistic regression beyond linear boundaries while maintaining proper uncertainty quantification.
Figure 15.1: Visualization of how neural networks approximate complex functions through composition of simple nonlinear transformations. The figure demonstrates a neural network approximating a target function (red solid line) composed of multiple sine waves of different frequencies. The network output (blue dashed line) closely matches the target pattern by combining multiple sigmoid neurons (light blue lines). Each individual neuron has a simple sigmoidal shape defined by a weight ww and bias bb as σ(wx+b)\sigma(wx+b) , but when linearly combined with learned coefficients, they form a powerful approximator capable of representing complex nonlinear relationships. This illustrates the principle behind neural networks: while individual neurons perform only simple nonlinear transformations, their composition creates models with extraordinary flexibility. Unlike polynomial regression which requires manual feature engineering, or Gaussian Processes which use kernel functions, neural networks automatically learn the appropriate feature transformations directly from data through the weights and biases of each neuron.
Figure 15.2: Evolution of activation functions in neural networks, highlighting the progression from early discrete models to modern differentiable alternatives. The Step function (red) was used in the original Perceptron (1958), creating a binary output that enabled basic linear classification but prevented gradient-based learning. The Sigmoid function (blue) and Tanh function (green), introduced in the 1980s, enabled backpropagation through their smooth, differentiable nature, though they suffered from vanishing gradient problems in deep networks. ReLU (purple), which became dominant in the 2010s, addressed this issue through its simple piecewise linear form that prevents saturation for positive inputs while providing true sparsity for negative inputs. GELU (orange), introduced in 2016, has become popular in modern transformer architectures, combining the benefits of ReLU with smooth transitions. This evolution reflects the requirement that introduced nonlinearity between linear transformations, which is essential for enabling neural networks to approximate complex functions through composition of simple operations. Each activation function represents a different trade-off between mathematical properties (differentiability, saturation, sparsity), with modern functions optimized for gradient-based learning in deep architectures.
Figure 15.3: Architecture of a fully connected feedforward neural network, also known as a Multi-Layer Perceptron (MLP), with input, hidden, and output layers. Each circle represents a neuron, with connections showing the information flow from inputs ( xix_{i} ) to outputs ( yjy_{j} ). Every neuron in one layer connects to every neuron in the subsequent layer, creating the “fully connected” structure. The magnified view shows the mathematical operation inside each neuron: a weighted sum of all inputs from the previous layer ( wTx+b\mathbf{w}^{T}\mathbf{x}+b ), passed through a nonlinear activation function ff . During the forward pass, information propagates from left to right as each layer transforms the outputs from the previous layer. This composition of simple transformations enables the network to learn complex representations—early layers detect basic features while deeper layers combine these features into higher-level abstractions. The key insight of neural networks lies in this hierarchical composition: by stacking multiple nonlinear transformations, the network can approximate virtually any continuous function to arbitrary precision, as guaranteed by the Universal Approximation Theorem.
Figure 15.4: Comparison between classical statistical learning theory and the modern double descent phenomenon observed in neural networks. The top panel shows the traditional U-shaped risk curve that has guided model selection for decades. In this classical view, as model complexity increases, test error initially decreases (addressing underfitting), reaches an optimal point, and then increases indefinitely as the model overfits. The bottom panel illustrates the double descent phenomenon, where test error follows a more complex pattern: it first decreases to a local minimum (similar to the classical curve), then increases as the model approaches the interpolation threshold (where parameters equal training samples), but—contrary to classical theory—begins decreasing again in the overparameterized regime, eventually achieving lower error than the first minimum. This unexpected asymptotic improvement in highly overparameterized models helps explain why modern neural networks with millions or billions of parameters can generalize well despite having far more parameters than training examples. This phenomenon has challenged conventional statistical wisdom and sparked new theoretical frameworks to understand the generalization capabilities of deep learning models.
Figure 15.5: Stochastic Gradient Descent (SGD) with different learning rates navigating a non-convex loss landscape. While this visualization was originally introduced in our discussion of logistic regression, the same principles apply directly to neural network training. The contour plots represent a simplified 2D projection of a high-dimensional loss landscape, with darker regions indicating lower loss values. The left panel shows SGD with a small learning rate, which follows the gradient more precisely but may become trapped in local minima. The right panel demonstrates how a larger learning rate enables the optimization to escape shallow local minima, potentially reaching better global solutions. In neural networks, this behavior is particularly important due to the highly non-convex nature of the loss landscapes. The stochasticity introduced by mini-batch sampling—as discussed in the text—provides an additional mechanism that works in concert with the learning rate to help navigate these complex loss surfaces.
Figure 15.6: Common learning rate schedules used in neural network training. The left panel shows inverse time decay, where the learning rate gradually decreases according to η(t)=η0/(1+γt)\eta^{(t)}=\eta_{0}/(1+\gamma t) . This schedule enables larger initial steps for faster progress in early training, followed by progressively smaller steps for fine-tuning. The right panel illustrates cosine annealing with restarts, where the learning rate oscillates according to η(t)=ηmin⁡+12(η0−ηmin)(1+cos⁡(tπ/T))\eta^{(t)}=\eta_{\min}+\frac{1}{2}(\eta_{0}-\eta_{\text{min}})(1+\cos(t\pi/T)) . This cyclical approach is particularly effective for neural networks, as it creates alternating phases of exploration (high learning rate) and exploitation (low learning rate), helping the optimization process escape local minima in complex loss landscapes. The restart mechanism, where the learning rate periodically returns to its maximum value, allows the optimizer to explore different regions of the parameter space.
Figure 15.7: Comparison of Stochastic Gradient Descent (SGD, left) and Adaptive Moment Estimation (Adam, right) on a loss landscape with varying curvature. Both optimizers start from the same point (black star) and attempt to find the global minimum (black circle). SGD uses a fixed learning rate that causes oscillatory behavior in regions of high curvature and prevents precise convergence to the minimum. In contrast, Adam combines momentum (shown by purple arrows) with adaptive learning rates to navigate the loss landscape more effectively. The momentum term m^t\hat{m}_{t} acts as a moving average of past gradients, helping Adam maintain consistent progress through regions with noise or high oscillation. Meanwhile, the second moment term v^t\sqrt{\hat{v}_{t}} scales each parameter’s update according to its gradient history, automatically taking smaller steps in directions with high gradient variance. Dividing the momentum term by the second moment estimate creates these adaptive learning rates, allowing Adam to make rapid progress in flat regions while precisely converging to minima without overshooting.
Figure 15.8: Visualization of gradient approximation using the finite difference method for a one-dimensional function f(x)=0.5x2−0.3sin⁡(3x)f(x)=0.5x^{2}-0.3\sin(3x) . At the point x0x_{0} , the true gradient (green dashed line) is given by the analytical derivative evaluated at that point. The finite difference approximations (red lines) estimate this gradient by computing f(x0+Δ)−f(x0)Δ\frac{f(x_{0}+\Delta)-f(x_{0})}{\Delta} for different values of Δ\Delta . As shown, the approximation becomes more accurate as Δ\Delta decreases, converging to the true gradient as Δ\Delta approaches zero. This illustrates the definition of derivatives: dfdx=lim⁡Δ→0f(x+Δ)−f(x)Δ\frac{df}{dx}=\lim_{\Delta\rightarrow 0}\frac{f(x+\Delta)-f(x)}{\Delta} . While conceptually simple and easy to implement, this approach becomes computationally prohibitive for neural networks with millions of parameters, as it requires a separate function evaluation for each parameter. For a network with PP parameters, at least P+1P+1 forward passes would be needed to compute the full gradient vector—one baseline evaluation plus one additional evaluation per parameter. This inefficiency motivates the development of backpropagation, which computes exact gradients more efficiently by applying the chain rule systematically through the computational graph.
Figure 15.9: Computational graph for a simple neural network showing the forward pass with stored intermediate values. The graph represents the function L=12(y−y^)2=12(y−σ(Wx+b))2L=\frac{1}{2}(y-\hat{y})^{2}=\frac{1}{2}(y-\sigma(Wx+b))^{2} decomposed into elementary operations. During the forward pass, computation flows left to right through the graph, with each node computing and storing its output value (shown in green boxes). These stored values—such as z=0.6z=0.6 at the linear computation stage and h=0.646h=0.646 at the sigmoid activation—are essential for efficient gradient calculation during backpropagation. This example illustrates the first phase of backpropagation, where the model establishes the specific point in parameter space at which gradients will be evaluated. The subsequent backward pass would propagate gradients from the loss node back to each parameter, applying the chain rule at each step to compute partial derivatives like ∂L∂W\frac{\partial L}{\partial W} and ∂L∂b\frac{\partial L}{\partial b} . By structuring computation as a graph with cached intermediate values, backpropagation achieves exact gradient computation with substantially greater efficiency regardless of the number of parameters and exact architecture of the neural network.
Figure 15.10: Visualization of the backward pass in backpropagation, showing how gradients flow from the loss function back to the network parameters. Following the forward pass (shown with the stored values in green), the backward pass begins at the loss node LL with ∂L∂L=1\frac{\partial L}{\partial L}=1 and propagates gradients right-to-left through the computational graph (red arrows). At each node, the incoming gradient is multiplied by the local derivative according to the chain rule. For instance, the gradient at y^\hat{y} is ∂L∂y^=−0.354\frac{\partial L}{\partial\hat{y}}=-0.354 , which then flows to hh with the same value since ∂y^∂h=1\frac{\partial\hat{y}}{\partial h}=1 . The gradient further propagates to zz as ∂L∂z=0.0697\frac{\partial L}{\partial z}=0.0697 after multiplication by ∂h∂z=σ′(z)=0.228\frac{\partial h}{\partial z}=\sigma^{\prime}(z)=0.228 . From zz , the gradient branches to compute ∂L∂W=0.0348\frac{\partial L}{\partial W}=0.0348 , ∂L∂b=0.0697\frac{\partial L}{\partial b}=0.0697 , and ∂L∂x=0.0697\frac{\partial L}{\partial x}=0.0697 . The detailed calculation for ∂L∂W\frac{\partial L}{\partial W} shown at the bottom illustrates how the full chain rule is applied: ∂L∂W=∂L∂y^⋅∂y^∂h⋅∂h∂z⋅∂z∂W\frac{\partial L}{\partial W}=\frac{\partial L}{\partial\hat{y}}\cdot\frac{\partial\hat{y}}{\partial h}\cdot\frac{\partial h}{\partial z}\cdot\frac{\partial z}{\partial W} . This systematic backward flow of gradients using stored forward pass values is what makes backpropagation computationally efficient, requiring only one forward and one backward pass regardless of the number of parameters and exact architecture of the neural network.
Figure 15.11: Comparison of forward mode and backward mode gradient calculation methods, illustrating why backpropagation is computationally efficient. The left panel shows forward mode calculation, where gradient dimensions grow exponentially with network depth, starting from inputs x∈RD\mathbf{x}\in\mathbb{R}^{D} and weights W∈RM×D\mathbf{W}\in\mathbb{R}^{M\times D} . The intermediate tensor ∂z∂W∈RM×(M×D)\frac{\partial\mathbf{z}}{\partial\mathbf{W}}\in\mathbb{R}^{M\times(M\times D)} rapidly expands in dimension. The right panel demonstrates backpropagation (backward mode), where gradient computation starts from the scalar loss L∈RL\in\mathbb{R} and propagates backward through ∂L∂y^\frac{\partial L}{\partial\hat{y}} , ∂L∂h\frac{\partial L}{\partial\mathbf{h}} , and ∂L∂z\frac{\partial L}{\partial\mathbf{z}} to ultimately compute ∂L∂W∈RM×D\frac{\partial L}{\partial\mathbf{W}}\in\mathbb{R}^{M\times D} . This approach maintains manageable dimensions throughout, resulting in reduced computational complexity.
Figure 15.12: Schematic comparison of Convolutional Neural Networks (CNNs) and Transformer architectures, highlighting their distinct structural designs and inductive biases. Top: The CNN processes an image through a sequence of specialized layers: convolution layers with sliding kernels detect local features, pooling layers reduce spatial dimensions while preserving important features, and fully connected layers integrate these features for the final output. This architecture encodes translation invariance and hierarchical local structure. Bottom: The Transformer processes sequential data through a different approach. Input tokens are embedded into a higher-dimensional space, then processed through a self-attention layer where each token interacts directly with every other token regardless of position (visualized by the curved connections). This global interaction is followed by a feed-forward network and output layer. The self-attention mechanism enables Transformers to capture long-range dependencies, making them particularly effective for data where relationships between distant elements are important, such as spectroscopic data where correlations between widely separated spectral lines can indicate specific physical conditions.
Figure 15.13: Comparison of sigmoid and ReLU activation functions and their derivatives. Top left: The sigmoid function σ(x)=1/(1+e−x)\sigma(x)=1/(1+e^{-x}) maps inputs to the range (0,1)(0,1) , but saturates for large positive or negative inputs (highlighted in red), causing vanishing gradients. Top right: The ReLU function ReLU(x)=max⁡(0,x)\text{ReLU}(x)=\max(0,x) and its smoother variant GELU (dashed purple) provide linearity for positive inputs while creating sparsity for negative inputs. Bottom left: The sigmoid derivative dσ(x)dx=σ(x)(1−σ(x))\frac{d\sigma(x)}{dx}=\sigma(x)(1-\sigma(x)) rapidly approaches zero in saturation regions, severely limiting gradient flow during backpropagation. Bottom right: The ReLU derivative is exactly 1 for all positive inputs, allowing strong and constant gradient flow, while being zero for negative inputs (which can lead to “dying neurons”). GELU provides a smooth transition around zero, addressing ReLU’s non-differentiability while maintaining its desirable properties for positive inputs. The vanishing gradient problem with sigmoid activations slows learning in deep networks, while GELU has emerged as a preferred choice in modern neural networks by combining ReLU’s advantages with improved smoothness and training stability.
Figure 15.14: Visualization of superlinear versus non-superlinear functions. The light red regions represent superlinear behavior where ∣y∣>∣x∣|y|>|x| , while the light blue regions represent non-superlinear behavior where ∣y∣≤∣x∣|y|\leq|x| . Five functions are displayed: linear ( f(x)=xf(x)=x ), quadratic ( f(x)=x2f(x)=x^{2} ), cubic ( f(x)=x3f(x)=x^{3} ), sigmoid, and GELU. Notably, quadratic and cubic functions quickly enter the superlinear regions, growing much faster than linear as ∣x∣|x| increases. In contrast, neural network activation functions (sigmoid and GELU) remain within the non-superlinear regions—they either saturate to a constant value (sigmoid) or grow at most linearly (GELU). Superlinear functions would amplify small input errors and cause unstable gradients in deep networks, making training difficult or impossible. The non-superlinear constraint represents a careful balance between expressivity and stability in neural network design.
Figure 15.15: Visualization of the dropout technique for neural networks. The left panel shows a standard neural network with fully connected layers, while the right panel demonstrates the same network with dropout applied (probability p=0.5p=0.5 ). During training with dropout, neurons are randomly deactivated (marked with red X’s) with probability pp , and their connections are temporarily removed. This forces the network to develop redundant pathways, preventing overfitting and creating an implicit ensemble of networks. Mathematically, dropout applies a random Bernoulli mask mi∼Bernoulli(1−p)m_{i}\sim\text{Bernoulli}(1-p) to each neuron and scales the outputs by 1/(1−p)1/(1-p) to maintain consistent expected values. When used during inference (Monte Carlo dropout), this technique approximates Bayesian neural networks by sampling from an implicit posterior distribution over network parameters, enabling uncertainty quantification in predictions.
Figure 15.16: Schematic illustration of the “grokking” phenomenon in neural networks. The top panel shows a conceptual representation of accuracy curves: while training accuracy (solid blue) quickly rises as the network memorizes examples, test accuracy (dashed red) remains low for an extended period before suddenly improving at the “grokking point.” The bottom panels illustrate the hypothesized transformation in neural representation: during memorization (left), neurons exhibit chaotic, seemingly random activation patterns, while after grokking (right), activations self-organize into structured patterns that reflect the underlying mathematical relationships in the data. This schematic visualization captures the essence of how neural networks can transition from memorizing specific examples to discovering generalizable principles. While simplified, this illustration represents patterns observed in research on neural networks trained on tasks like modular arithmetic, where systematic analysis of internal representations reveals that networks eventually discover compact, generalizable solutions rather than merely interpolating between training examples. This phenomenon challenges the view of neural networks as fancy interpolators and suggests they can develop representations that capture the abstract structure governing the data.
Figure 15.17: Schematic illustration comparing foundation models with few-shot learning to traditional training approaches. The top panel presents a conceptual performance curve showing how foundation models (blue) achieve high accuracy with very few examples, while models trained from scratch (red) require orders of magnitude more data to reach comparable performance. The shaded blue area highlights the performance advantage of foundation models, particularly pronounced in the few-shot region. The plot demonstrates how performance varies across few-shot, small data, and large data regimes. The bottom panel diagrams the foundation model paradigm: a foundation model is pre-trained once on diverse, extensive data (left blue box), then adapted to specific tasks with minimal examples through fine-tuning (top right green box), resulting in specialized models (right orange box). This approach contrasts with the traditional method of training separate models from scratch for each task (shown in gray). This conceptual diagram illustrates the computational efficiency of the foundation model approach, where the high computational cost of pre-training is incurred only once, while adaptation requires minimal computation and can be performed quickly for new tasks.
Figure 15.18: Architecture of an autoencoder neural network for nonlinear dimension reduction. The hourglass structure illustrates the compression and reconstruction process: the encoder (blue) progressively reduces input dimensionality from 1000 to 16 through hidden layers of decreasing size (256, 64), while the decoder (green) symmetrically expands from the bottleneck layer back to the original dimensionality. The bottleneck layer (orange) creates an information bottleneck that forces the network to learn efficient representations of the data. During training, the network is optimized to minimize reconstruction error between the input x\mathbf{x} and output x^\hat{\mathbf{x}} , encouraging the latent space z\mathbf{z} to capture meaningful features. Unlike Principal Component Analysis (PCA) which is limited to linear transformations, autoencoders can model complex nonlinear relationships through their neural network layers.
Figure 15.19: Comparison of linear and nonlinear dimensionality reduction techniques using the Swiss roll dataset. The left panel shows the original Swiss roll data—a two-dimensional manifold (sheet) curled in three-dimensional space, with colors indicating positions along the spiral. The middle panel demonstrates Principal Component Analysis (PCA), a linear technique that projects data onto directions of maximum variance. PCA fails to properly unfold the Swiss roll structure because it is constrained to linear transformations, resulting in overlapping points from different parts of the manifold (note the mixing of colors). The right panel illustrates the result of an autoencoder-based dimensionality reduction. By leveraging nonlinear transformations, the autoencoder effectively “unrolls” the Swiss roll, preserving local neighborhood relationships and maintaining the continuous color progression that represents positions along the original manifold.
Figure 15.20: Encoder-decoder architectures for transformations between data domains. (a) Domain transfer maps data between different representations of the same type, such as converting simulated stellar spectra to match observed spectra. This approach helps bridge the gap between theoretical models and observations by learning systematic differences between simulations and real data. (b) Shared latent space architectures enable integration of different data modalities by mapping them to a common representation. For example, by encoding both galaxy images and spectra into the same latent space, the network learns correlations between visual morphology and spectral features, forcing the network to discover relationships between different data types.
Figure 15.21: Comparing deterministic and probabilistic approaches to modeling one-to-many relationships. Panel (a) shows how a standard neural network trained with mean squared error (red curve) models the Swiss roll dataset. This deterministic approach must predict a single output for each input, effectively averaging across multiple valid solutions at the same xx coordinate. Panel (b) shows a probabilistic approach with a vertical dashed line at x=2.5x=2.5 demonstrating how one input can correspond to multiple possible outputs. The right panel shows the conditional distribution p(y∣x=2.5)p(y|x=2.5) , with multiple peaks capturing all possible outcomes and their probabilities. Mixture Density Networks model the full conditional distribution p(y∣x)p(y|x) that varies smoothly with xx , allowing the model to capture how the multimodal nature of the data gradually changes across the input space.
Figure 15.22: Architecture of a Mixture Density Network (MDN) for modeling conditional probability distributions. The left side shows the neural network structure with input neurons (blue), hidden neurons (green), and output neurons organized by mixture component (orange, purple, red). Each mode has three parameter outputs: mixture weight ( πk\pi_{k} ), mean ( μk\mu_{k} ), and variance ( σk2\sigma^{2}_{k} ). The right side visualizes the resulting conditional distribution p(y∣x)p(y|\mathbf{x}) as a mixture of Gaussians, with each colored component corresponding to a mode in the network output. MDNs extend traditional neural networks to model full probability distributions rather than point estimates, making them particularly valuable for problems with one-to-many relationships where a single input can correspond to multiple plausible outcomes. The network learns to output the parameters of a mixture model that varies with the input, effectively implementing the conditional distribution p(y∣x)=∑k=1Kπk(x)N(y∣μk(x),σk2(x))p(y|\mathbf{x})=\sum_{k=1}^{K}\pi_{k}(\mathbf{x})\mathcal{N}(y|\mu_{k}(\mathbf{x}),\sigma^{2}_{k}(\mathbf{x})) . This approach enables uncertainty quantification and multimodal predictions, addressing a limitation of standard (deterministic) neural networks in scientific applications.
Figure 15.23: Visualization of normalizing flows showing the gradual transformation between probability distributions. The diagram illustrates how a complex bimodal target distribution (left, red) is progressively transformed through intermediate steps into a standard Gaussian base distribution (right, blue). The gray lines track specific points ( x1x_{1} through x5x_{5} ) as they move through each transformation stage, demonstrating how probability mass is transported from one space to another. The forward direction (left to right) enables density evaluation using the change of variables formula shown at the bottom, while the inverse direction (right to left) allows sampling from the complex target distribution. The Jacobian determinant in the formula accounts for how volume elements are stretched or compressed during transformation. Normalizing flows learn this invertible mapping directly from data, providing both precise density estimation and efficient sampling capabilities—a combination that makes them particularly powerful for modeling complex astronomical distributions where analytical forms for likelihoods are often intractable.
Figure 15.24: Visualization of diffusion models showing the gradual addition of noise to transform data into a standard Gaussian distribution. The diagram illustrates how a complex bimodal data distribution (left, dark blue) progressively absorbs more noise through intermediate time steps ( t=0.2,0.5,0.8t=0.2,0.5,0.8 ) until reaching pure Gaussian noise (right, light blue). Specific points ( x1x_{1} through x5x_{5} ) are tracked throughout the process, eventually becoming xT,1x_{T,1} through xT,5x_{T,5} at the pure noise stage. Small 1D Gaussian distributions alongside each point visually represent how noise progressively disperses the data, with the width of these mini-Gaussians increasing with each step to show the growing variance. The forward process (dark blue arrow) follows a Markov chain of noising steps defined by the equation xt=1−βtxt−1+βtεx_{t}=\sqrt{1-\beta_{t}}x_{t-1}+\sqrt{\beta_{t}}\varepsilon , where βt\beta_{t} controls the noise schedule. The learned reverse process (light blue arrow) performs denoising to generate new samples. Unlike normalizing flows which require invertible transformations, diffusion models operate through this progressive noising-denoising approach, allowing them to model complex distributions without architectural constraints. This approach has proven particularly effective for generating high-quality astronomical data where physical consistency and uncertainty quantification are essential.
Figure 15.25: Visualization of Simulation-Based Inference (SBI) showing the relationship between parameter space and observation space. The central panel shows the joint distribution p(θ,x)p(\theta,x) with black dots representing individual simulations (parameter-observation pairs). The top panel displays Neural Posterior Estimation (NPE), showing the posterior distribution p(θ∣x)p(\theta|x) obtained by conditioning on a specific observation value xx (purple horizontal dashed line). The right panel shows Neural Likelihood Estimation (NLE), displaying the likelihood function p(x∣θ)p(x|\theta) obtained by conditioning on a specific parameter value θ\theta (red vertical dashed line). Both approaches enable amortized inference—once trained, they can rapidly evaluate likelihoods or posteriors for any values without re-running simulations. This visualization highlights how SBI addresses a challenge in astronomical inference: NPE directly learns posteriors in lower-dimensional parameter space (more computationally efficient), while NLE maintains the traditional likelihood-based framework but works with complex simulation-based models. By leveraging neural networks to learn these relationships from simulations, SBI enables Bayesian inference for sophisticated models where analytical likelihoods would be intractable, without sacrificing statistical rigor.
More in this topic

Supervised Learning

See all 31 papers

Explore similar work

CardsList
Sep 18, 2026

Bayesian classification of astronomical spectra with class uncertainties

Context: We developed a probabilistic machine learning method with the aim of performing the O(10)-way classification of low- and high-resolution spectra of stellar and extragalactic targets for the upcoming 4MOST survey. In fulfilment of the survey requirements, this method should be able to express uncertainty in the input data as well as uncertainty introduced in its prediction. Aims: Four different methods are explored: (1) convolutional neural networks (CNNs), (2) the Dirichlet distribution, (3) Monte Carlo dropout (MCD), (4) Bayesian neural Networks (BNNs) + variational inference (VI). Training and validation was performed using labelled spectra from the SDSS database and a custom 4MOST mock dataset. All the methods were compared in terms of the same metrics: accuracy, area under the curve (AUC), expected calibration error (ECE), Shannon entropy, negative log-likelihood (NLL), Brier score, training time, and inference time. Methods: A CNN with simple architecture and about 20,000 parameters was trained to achieve classification accuracies of 91.5% on SDSS data and 92.8% on 4MOST mock data. The direct Dirichlet prediction and VI models tested provide uncertainties on class membership probabilities, but they confuse classes more often. The MCD on a CNN is found to be the most suitable; it boosts the point-estimate accuracies to 92.6% and 93.9%, while still providing fast training and sufficiently fast inference. Compared to a standard CNN, the method additionally provides well-calibrated uncertainties at marginal extra cost.

Simon Barton, Martin Sahlén, Andreas Korn +1

Theoretical Astrophysics, Division of Astronomy and Space Physics, Department of Physics and Astronomy, Uppsala University, Box 524, 751 20 Uppsala, Sweden. · Observational Astrophysics, Division of Astronomy and Space Physics, Department of Physics and Astronomy, Uppsala University, Box 524, 751 20 Uppsala, Sweden · Division of High Energy Physics, Department of Physics and Astronomy, Uppsala University, Box 524, 751 20 Uppsala, Sweden +1

Monte Carlo DropoutMulticlass ClassificationUncertainty CalibrationBayesian Neural Networks
AskPDF
May 18, 2026

Hyrax: An Extensible Framework for Rapid ML Experimentation and Unsupervised Discovery in the Era of Rubin, Roman, and Euclid

The NSF-DOE Vera C. Rubin Observatory, Roman Space Telescope, Euclid, and other next-generation surveys will deliver imaging, spectroscopic, and time-domain data at scales that increasingly shift the bottleneck in astronomical machine learning (ML) projects from model design to infrastructure. We present Hyrax, an open-source, modular, GPU-enabled Python framework that supports the full ML lifecycle in astronomy: from data acquisition and training to inference and experiment comparison, with capabilities including multimodal dataset support, integrated vector databases for similarity search, and interactive two- and three-dimensional latent-space exploration for unsupervised discovery. We demonstrate Hyrax's versatility through five representative applications on real survey data: (i) unsupervised representation learning on ∼4×105\sim 4\times10^5 Rubin Legacy Survey of Space and Time (LSST) Data Preview 1 (DP1) galaxies, surfacing new merger and low-surface-brightness candidates missing from reference Euclid and Dark Energy Survey catalogs, while also isolating imaging artifacts -- all without labeled training data; (ii) hybrid density-based clustering for identifying cluster-scale gravitational lens candidates in DP1 data; (iii) multimodal early-time transient classification in the Zwicky Transient Facility leveraging light curves, spectra, images, and metadata; (iv) supervised false-positive filtering in shift-and-stack searches for distant solar system objects in the Dark Energy Camera Ecliptic Exploration Project survey; and (v) supervised detection of semi-resolved dwarf galaxies in Hyper Suprime-Cam and LSST-like imaging using synthetic source injection. Together, these results demonstrate that Hyrax provides astronomy-specific ML infrastructure that enables systematic discovery and rapid methodological iteration across next-generation astronomical surveys.

Aritra Ghosh, Drew Oldag, Michael Tauraso +29

Dept. of Astronomy & the DiRAC Institute, University of Washington, Box 351580, Seattle, WA 98195, USA · eScience Institute, University of Washington, Seattle, WA, U.S.A. · School of Physics and Astronomy, University of Minnesota, Minneapolis, MN 55455, USA +6

Scientific MLAstrophysics
AskPDF
Jun 9, 2026

Interpretable Neural Marked Statistics for Cosmological Inference

Recovering cosmological information beyond the power spectrum is a central goal for upcoming cosmological surveys, since late-time non-Gaussian signal in the matter density cannot be accessed through two-point statistics alone. Marked statistics fold part of this information back into the two-point level by reweighting the field with non-linear functions. We propose a neural marking scheme to generalize this process through a set of interpretable, physically motivated transformations that directly allow to interpret the gain in cosmological information at the morphological level. We employ a contrastive learning objective to align learnable marked summaries with the underlying cosmological parameters. At kmax⁡=0.2 hMpc−1k_{\max}=0.2\,h\mathrm{Mpc}^{-1}, our neural mark tightens the marginalized constraint on σ8σ_8 by 2.9×2.9\times and on ΩmΩ_m by 1.8×1.8\times compared to classical marks, breaking the Ωm−σ8Ω_m-σ_8 degeneracy at the Fisher information level. It further reduces the parameter MSE across our cosmological parameter prior by 1.45×1.45\times over the best classical mark. The learned latent geometry aligns with the ΩmΩ_m and σ8σ_8 directions in parameter space, indicating that the contrastive objective recovers the dominant axes of cosmological information. Our approach opens the door to more powerful, interpretable summary statistics for cosmological inference.

Federico Semenzato, Benjamin D. Wandelt, Michele Liguori +1

Dipartimento di Fisica Galileo Galilei, Universit`a di Padova, I-35131 Padova, Italy · INFN Sezione di Padova, I-35131 Padova, Italy · Department of Physics and Astronomy, Johns Hopkins University, 3400 North Charles Street, Baltimore, MD, 21218, USA +2

Cosmology
AskPDF