Identifying Undiscovered Millisecond and Young Pulsars with Bayesian Neural Networks
Advisor: Bhattacharya
Internship Location: San Francisco State University
Onsite Mentor: Dr. Oscar Macías, Assistant Professor of Multi-Messenger Astronomy, San Francisco State University
Pulsars are rapidly rotating, highly magnetized neutron stars that emit beams of electromagnetic radiation, observed as periodic pulses when those beams reach Earth. In modern radio surveys, the candidate selection problem arises because millions of pulsar-like signals are detected, making it challenging to reliably identify true pulsars. My project focuses on solving the candidate selection problem of pulsars through the use of a Bayesian neural networks (BNN).
With the help of the Department of Physics & Astronomy at San Francisco State University, spectral energy distribution data will be gathered from the most recent Fermi Large-Area Telescope and pulsar catalogs. These will be used as the BNN's input data to classify millisecond and young pulsars with reliable probability distribution estimates. Based on these uncertainty quantifications, other researchers will be able to determine which high-confidence candidates to pursue in future studies. This will be significant in discovering new objects in our universe and understanding the Galactic structure.
Project Blog
Week 1: Prepare for Liftoff
Hi everyone! Welcome to week one of my senior project blog posts! My name is Chloe, and I’m super excited to take you along with me through all of the feats and flops of my project over the next several weeks. This week, I'll start light and easy so as not to scare any of my readers away. Today, I’ll introduce my project and optimistically describe my goals. So please, grab yourself a cup of coffee or tea, get comfortable, and prepare for liftoff.
For the next several weeks, I’ll use machine learning algorithms to automate and optimize what astronomers call the “pulsar candidate selection” process. These candidates refer to a vast array of possible pulsars, which are highly magnetized, rapidly rotating neutron stars formed from supernova explosions. They emit periodic pulses of radiation, which is how they get their name, and these intervals can range anywhere between seconds and milliseconds!
These cosmic objects are exotic and fascinating to study, but they are also extremely difficult to identify and categorize, due to the noisy nature of massive radio telescope datasets. As modern telescopes produce millions of signals, manual detection is nearly impossible (I’m shuddering just thinking about sitting in front of ten laptops and straining my eyes on hundreds of noise spikes). Researchers have been developing automated machine learning models to distinguish real pulsars from background noise. With these new optimized tools, a process that once took days can now be reduced to mere seconds.
So...What's the Point?
My project rides on this new wave of automated methods, focusing specifically on the classification between pulsars and other gamma-ray sources in the center of the Milky Way Galaxy. By gamma-ray sources, I mean the most energetic objects in the universe. These objects produce radiation with the smallest wavelengths, often involving high-energy processes. Some examples include black holes, neutron stars, supernova remnants, and pulsars, which are the main focus of my project. Differentiating between pulsars and other types of gamma-ray sources is crucial to understanding the structure of the universe.
With so much talk about pulsars, you might be wondering, “What different types of pulsars even are there?” Fear not, my attentive reader, for I am here to relieve all of your confusion. Today, I will be focusing on the two main types of pulsars I aim to classify: young and millisecond pulsars. While both types of pulsars are extremely exotic and energetic objects, they differ in age, spin speed, and magnetic field strength. These differences are necessary to understand the origins of magnetic fields and pulsar activities.
Young pulsars are, quite literally, young. They are newly formed pulsars often found in the center of the Galaxy, near supernova remnants. They have strong magnetic fields and a large amount of energy, but they also lose this energy at a much faster rate. On the other hand, millisecond pulsars are much older, “recycled” pulsars whose periods are in the millisecond range. They start as young pulsars, but are later spun up by accreting material from a companion star. Compared to young pulsars, millisecond pulsars have much weaker magnetic fields and slower rates of rotation decay, making them much more stable.
Looking Forward
Now that you've hopefully gotten a better understanding of pulsars, I can introduce you to my project's ultimate goal. My project aims to contribute to the ongoing endeavor of automated methods by using a Bayesian neural network. Unlike traditional neural networks that provide a simple binary classification, Bayesian neural networks provide an additional layer of information through probability distributions, allowing us to quantify the model’s uncertainty. So, not only will the model reveal what type of gamma-ray source we’re studying, but it will also show how confident it is in its own predictions!
Using these confidence measures, I will create a ranked list of pulsar candidates for researchers to pursue in future studies. With high-confidence pulsars at the top, this list will hopefully act as a foundational benchmark in pulsar research. Just like Neil Armstrong said, even if it’s just a small step for one high school senior, it could be a big leap for pulsar discovery.
This week, I’ve officially started my journey by gathering and organizing my data. Given my position as a high school senior with a limited budget, I won’t be using any fancy, expensive instruments to launch my own space telescope. Instead, I’m using two of NASA’s public Fermi Large Area Telescope catalogs–a comprehensive 14-year catalog[1] with many types of energetic, astronomical objects, and another catalog focusing on pulsars[2] –that can be found online to scrape all of my numerical data. After doing some data mining and preprocessing, I’ve been able to curate a clean dataframe of over 200 pulsars.
In the next couple of weeks, I’ll be using this data to develop a reliable framework for classifying pulsars. This could advance our understanding of high-energy astrophysical phenomena and contribute to machine learning in modern astronomy. I’m sure there will be so many ups and downs throughout this project, but I’m excited to show you the unfiltered journey. And with that, I must bid you farewell for now. I hope to see you very, very soon in next week’s blog post!
References
[1] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
[2] NASA Fermi-LAT Collaboration. (2025). LAT Third Catalog of Gamma-ray Pulsars. Version 3, Fermi Science Support Center, 21 Jul. 2025, fermi.gsfc.nasa.gov/ssc/data/access/lat/3rd_PSR_catalog/. [Dataset].
[3] Berry, Dana. "Accretion of Material Spins Pulsar to Millisecond Range." Wikimedia Commons, 7 Sep. 2013, commons.wikimedia.org/wiki/File:Accretion_of_Material_Spins_Pulsar_to_Millisecond_Range.jpg. Public Domain.
[4] NASA Chandra Team. "A new Chandra movie of the Vela pulsar shows it may be “precessing,” or wobbling as it spins." NASA, 7 Jan. 2013, www.nasa.gov/image-detail/a-new-chandra-movie-of-the-vela-pulsar-shows-it-may-be-precessing-or-wobbling-as-it-spins/. Public Domain.
Week 2: Blast from the Past Literature
Hello, hello! Thanks for stopping by again to follow me through week two of my senior project. If you’re reading this, then I guess that means you’re ready to learn more about the progress I’ve made since last week’s blog post. Since you’ve so kindly requested it, I’m happy to report that I’ve got plenty of exciting progress and updates to share!
This week, I will focus primarily on describing my dataset and doing a mini literature review to paint the competitive landscape for my project. Yes…that means things might get a bit technical today, but don’t fret! Your faithful narrator is here to break it all down for you, step by step.
Please Welcome...the Dataset
To start, let’s continue last week’s discussion of my dataset. I’d like to explain in more detail the specific variables and pulsar characteristics I’ve collected. Using two of NASA’s Fermi Large Area Telescope catalogs[1,2], I compiled a dataset of exactly 294 pulsars, including both millisecond and young pulsars. Of course, the model needs more information than just the two pulsar names if it’s going to learn any meaningful patterns for classification. Thanks to the Fermi collaboration’s extensive data collection and calculations, I was able to mine 30 additional features from the catalogs. All of these features are numerical, and together they can capture each pulsar's identity and spectral profile.
By spectral profile, I mean the characteristics that describe a pulsar’s energy distribution across a range of energies. Each pulsar’s energy spectrum is essentially a snapshot of how much radiation it emits across an energy range. Astronomers usually refer to this spectrum as a spectral energy distribution (SED). Rather than looking at the object, a pulsar in our case, in a single energy band, the SED is essentially a curve showing how the pulsar's brightness varies across the entire energy spectrum. This curve can then be summarized using parameters like the spectral slope and where the emission begins to drop at higher energies. These parameters are important because they can reveal the physical processes within a pulsar that result in those specific spectrum curves. Because young and millisecond pulsars tend to produce subtly different spectral shapes, these spectral parameters can be a useful, direct piece of data to distinguish between the two classes.
Paper Trails
The great thing about using all of this numerical data is that it allows us to extract relevant information for the classification directly from the energy spectra. This means that we can avoid potentially losing any valuable information, which would have been inevitable had we instead derived measurements.
This decision is crucial in ensuring that our data is as robust and informative as possible. Previously, other researchers have tackled the pulsar classification problem with intricate and complex models. Eatough et al.[3] were the first to introduce the idea of using an artificial neural network (ANN) to classify pulsars using their candidate information. An ANN is essentially the most basic form of a neural network, designed to mimic the human brain’s structure by using layers of interconnected “neurons” to analyze data and recognize complex patterns. It learns by adjusting connections between these nodes based on the input data, so the choice of data is very important.
Much of the past research has relied heavily on input that was either derived or extracted from noisy signals[2,3,4]. Multiple papers have used signals detected from pulsars’ periodic pulses to train their machine learning models. While pretty effective, the approach requires working with observational data that can contain a lot of noise, which must be filtered before analysis.
Other work, such as that of Zhu et al.[5], has tried applying image-recognition techniques to identify pulsars based on visual plots of these pulse characteristics, produced during pulsar searches. This is a novel approach, but it still requires the model to learn from complex representations and noisy observations.
Marking My Territory
Unlike much prior research, my project simplifies this process by focusing directly on the numerical spectral parameters describing each pulsar’s gamma-ray energy spectrum. This will allow the model to learn from clean, structured data rather than from noisy pulses or visual plots. It will also be a lot easier to read the results and to identify notable pulsar differences, as we’ll be able to pinpoint exactly which range of numbers marks the strongest division between the two types of pulsars. Just think about how much more difficult it would be to interpret a model trained on complex candidate images!
Perhaps the biggest difference between my model and earlier models is its ability to quantify its own confidence. Whether it’s improving the neural network approach or introducing other types of machine learning models[6,7], much of the previous literature has relied solely on a single predicted class probability. While it’s already incredibly useful to have this as a final classification, adding a distribution of predicted probabilities lets us identify the pulsars our model classified with the highest confidence.
That is where the "Bayesian" part in my project’s title finally comes in. I’ll be discussing that important concept and explaining how it allows us to produce a whole range of predictions instead of just a single number. Until then, stay tuned!
References
[1] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
[2] NASA Fermi-LAT Collaboration. (2025). LAT Third Catalog of Gamma-ray Pulsars. Version 3, Fermi Science Support Center, 21 Jul. 2025, fermi.gsfc.nasa.gov/ssc/data/access/lat/3rd_PSR_catalog/. [Dataset].
[3] Eatough et al. “Selection of radio pulsar candidates using artificial neural networks.” Monthly Notices of the Royal Astronomical Society, vol. 407, no. 4, 2010, pp. 2443-2450. Oxford Academic, doi.org/10.1111/j.1365-2966.2010.17082.x.
[4] Bates et al. “The High Time Resolution Universe Pulsar Survey — VI. An artificial neural network and timing of 75 pulsars.” Monthly Notices of the Royal Astronomical Society, vol. 427, no. 2, 2012, pp. 1052–1065. Oxford Academic, doi.org/10.1111/j.1365-2966.2012.22042.x.
[5] Zhu et al. “Searching for Pulsars Using Image Pattern Recognition.” The Astrophysical Journal, vol. 781, no. 2, 2014, pp. 117-129. IOPscience, 10.1088/0004-637X/781/2/117.
[6] Devine et al. “Detection of dispersed radio pulses: a machine learning approach to candidate identification and classification.” Monthly Notices of the Royal Astronomical Society, vol. 459, no. 2, 2016, pp. 1519–1532. Oxford Academic, doi.org/10.1093/mnras/stw655.
[7] Lyon et al. “Fifty years of pulsar candidate selection: from simple filters to a new principled real-time classification approach.” Monthly Notices of the Royal Astronomical Society, vol. 459, no. 1, 2016, pp. 1104–1123. Oxford Academic, doi.org/10.1093/mnras/stw656.
Week 3: A Bayesian Test Flight
Hi everyone, and welcome back! This week, I'll finally get into the details of the “Bayesian” part of my project’s title. For those of you with some experience in statistics, this term might sound familiar. Indeed, Bayesian neural networks (BNNs) draw on the fundamental field of Bayesian statistics. Bayesian statistics is an approach to statistical inference grounded in Bayes’ theorem that updates the probability of a hypothesis as new evidence becomes available. Don’t worry if all of that sounded like academic word soup. This entire post will be dedicated to explaining the concept step by step in a much more approachable way!
I spent the entirety of this week poring over the contents of two papers: Jospin et al.’s “Hands-on Bayesian Neural Networks – A Tutorial for Deep Learning Users,”[1] which provides a thorough tutorial on the theory and algorithms behind BNNs; and Butter et al.’s “Classification of Fermi-LAT blazars with Bayesian neural networks,”[2] which puts the idea of a BNN into action. After having read through 51 pages of dense neural network sorcery, I’ve gained a solid understanding that I’ll now be able to convey to you, my brave and curious reader.
Reproducing Results
I won't simply regurgitate the equations and mathematical jargon I encountered. Instead, I’ll explain the Bayesian concepts more actively by recounting how I built my own BNN. Thanks to Butter and co’s well-explained work and methodology in their paper (which I highly encourage everyone to read!), I was able to create a similar BNN and reproduce some of their results.
I chose to begin by “copying” the results of Butter et al. because: 1) their research is quite similar to mine in that they employ BNNs to classify astronomical objects, and 2) one must learn the ropes before setting sail.
Naturally, I began by gathering the data used in Butter et al.’s paper. Although they were classifying blazars (no need to worry about what blazars are—just think of them as another class of energetic objects powered by supermassive black holes in space), they used the same Fermi Large Area Telescope (Fermi-LAT) catalog[3] I used to mine my pulsar data in last week’s discussion. Given my familiarity with the structure of the Fermi-LAT catalog, this step was easy enough.
Within a couple of minutes, I had gathered the same data as the paper. This data consisted of 2310 blazar sources, along with the energy spectrum of each blazar. The energy spectrum was reported as a series of fluxes, or intensities of radiation, across seven consecutive energy bands. These bands came from dividing the energy range 50 megaelectronvolts (MeV) to 300 gigaelectronvolts (GeV) into seven subranges (1: 50–100 MeV; 2: 100–300 MeV; 3: 300 MeV–1 GeV; 4: 1–3 GeV; 5: 3–10 GeV; 6: 10–30 GeV; 7: 30–300 GeV).
Following the paper’s methodology, I took the logarithms of these flux values. I then computed the median value for each energy band separately. The plot below shows the median energy spectra for the two blazar classes, BLL and FSRQ (again, just understand that these are two types of blazars). The connected lines represent the median log flux, and the shaded bands represent the lower and upper quartiles.
It Gets Math-y
With the dataset all ready to go, we’re ready to train the neural networks…right? Not so quick! We first need to understand neural networks and Bayesian statistics before we can confidently execute the methodology ourselves.
Let’s start with traditional dense neural networks (DNNs). These are also known as “fully connected layers,” which makes it a little easier to envision what they are. Just like their alternative name suggests, DNNs are a type of neural network where each neuron in one layer is connected to every neuron in the next layer. This fully connected structure of DNNs allows them to learn complex relationships in the data. The graphical representation I created below shows a simple DNN setup. Notice how each neuron (represented as a shaded circle) in each layer has a line going to every neuron in the following layer.
The first layer (on the left) is the input layer, which holds the raw data features the network must learn from. Right after the input layer come the hidden layers. This is where all the fun happens! Here, the data passes through these layers and is transformed by applying weights and biases that are adjusted during training. The weights are essentially just matrices that connect the neurons and determine the strength and direction of their connections, while the biases are constants added to shift the neuron’s input.
These weights and biases are the parameters the network must learn throughout the training process. In a dense neural network, the weights are simply single values within a matrix. With every run through the network, known as a forward pass, the neurons receive input values from each layer, multiply them by their weights, add the biases, apply activation functions to introduce non-linearity, and pass the outputs to the next layer. This is how the weights and biases get updated with every pass.
However, it wouldn’t be much use for the network to run once and call it a day, would it? It wouldn’t be possible to actually learn anything that way. That’s why networks use a handy optimization system called backpropagation and gradient descent to gradually adjust their internal parameters and improve their predictions over time. After every pass, the network calculates the difference between predicted and actual values to compute the loss. Backpropagation is then the process of working backwards through the network to determine how much each connection and parameter contributed to that error. Once this information is known, the network uses gradient descent to make tiny adjustments to those parameters in the direction that reduces the loss. The model repeats this cycle many, many times to improve its performance and make more accurate predictions.
This is essentially the basic training process for all types of neural networks. Different network variations add their own minor tweaks, but the overall idea is very similar. Now that we’ve all got a general idea of how neural networks work, we can return to recreating Butter et al.’s results. The authors decided that they would experiment with creating a DNN before jumping directly into building their final BNN. This baseline DNN would then be able to serve as a point of comparison for the BNN’s performance. A very sensible move, if you ask me! Following in their footsteps, I implemented the same DNN.
The authors very conveniently provided the exact architecture of their DNN: they used the seven energy bins with their corresponding log flux values as the input vector, one hidden layer, and 16 neurons in that hidden layer. I was easily able to recreate their exact DNN model, train it, and evaluate it using a standard 10-fold cross-validation procedure. To quantify the DNN’s performance, I created a receiver operating characteristic (ROC) curve, which plots the true positive rate against the false positive rate to show the performance across all possible classification thresholds. The plot is shown below with its area under curve (AUC), along with a confusion matrix of classified objects.
Predicted BLL | Predicted FSRQ
Actual BLL | 1375 | 115
Actual FSRQ | 143 | 677
Left: ROC curve of my recreated DNN with corresponding AUC value. Right: confusion matrix of classifications
These metrics served as a kind of starting point and allowed me to familiarize myself with constructing a DNN before tackling the more complex sibling: the BNN.
Bayesian Statistics
The main difference between a traditional neural network and a BNN is the structure of their weights. What I mean is that a simple DNN has single, point values for weights, while a BNN has distributions for weights. So rather than learning a single number for each weight, a BNN must learn the distribution of each weight, usually characterized by the weight’s mean (approximate value) and variance (spread or uncertainty in the prediction).
Training these weights relies on a fundamental idea from statistics known as Bayes’s theorem, which describes how we update beliefs as we see new evidence. The theorem is
P(θ | D) = P(D | θ) · P(θ)P(D)
Allow me to break down this seemingly intimidating equation. Each term in this equation has an intuitive meaning. Here, θ represents the hypothesis for the neural network’s weights, and D represents the training data. The quantity P(θ | D), called the posterior, is what we ultimately want the BNN to learn: it’s the updated distribution of weight values after the network has seen the data. The term P(D | θ) is called the likelihood and measures how well a particular set of weights explains the data. The factor P(θ) is called the prior, or what we initially assume about the weights prior to training. For example, we can often assume that the most reasonable weight values are centered near zero. The final term P(D) is just a normalization factor that makes sure the probabilities add up correctly.
In practice, training a BNN means finding the true posterior distribution over the weights that balances what the data suggests and what we assumed beforehand. However, directly calculating the true posterior is usually complex and computationally infeasible, so researchers use approximation techniques such as variational inference. Rather than computing the exact posterior, variational inference assumes a simpler family of variational distributions and then adjusts its parameters so that this approximation is as close as possible to the true posterior. Butter et al.’s paper makes use of this approach to efficiently train their BNN, so of course, I had to follow.
Battle of DNN vs. BNN
Using the same architecture as that of the DNN, I set out to recreate the BNN. The process was very similar, except this time I used probabilistic network layers to capture the weight distributions. I constructed the BNN and experimented with a toy dataset to compare the predictions of the BNN with those of an ensemble of independently trained DNNs. As Butter et al.’s paper discusses, both of these approaches would be able to estimate the uncertainty of a prediction, and it would be a good strategy to test the performance of each neural network.
Rather than testing both neural networks on the entire dataset, I followed Butter’s team’s instructions to create a toy dataset with a structure similar to the original dataset. I trained the BNN by sampling 100 different possible weight distribution configurations. I then trained an ensemble of 100 DNNs using two different methods: one where all 100 DNNs were trained on the same dataset, and another where 100 training datasets were sampled independently. This was a choice made by Butter and co to estimate the impact of reusing small training datasets.
After training each of the three models, I created plots to visualize the correspondence between their predictions. The plots below show the correlations between the resulting means and standard deviations of the BNN and DNN predictions.
From these plots, we observe a relatively linear correlation in the mean predictions. There also seems to be a stronger correlation between the BNN and the DNN trained on independent data than between the two DNNs. This implies that their performance was much more similar. The lower panels show the standard deviation correlations. We notice more interesting patterns over here: the DNN trained on independent samples significantly underestimates uncertainty compared to the BNN. Meanwhile, training multiple DNNs on the same dataset seems to underestimate the uncertainties. These plots show the need for BNNs to obtain reliable uncertainty estimates.
This entire experiment has solidified my familiarity and confidence in working with various neural network architectures. While I feel much more prepared, I must do much more model fine-tuning. For one, I must continue experimenting with the probabilistic layers to see if I can get predictions closer to those in Butter et al.'s paper. I must also become more experienced in interpreting the unique loss functions of BNNs.
I will document all of these obstacles and discoveries in my next blog post. But for now, I must bid you farewell again. Thank you for bearing with me through this admittedly very technical post. Now that we've got the complicated background knowledge out of the way, future posts will be much easier to digest. With that, I'll see you very, very soon!
Acronyms List
For your convenience, here is a list of all the acronyms. I ordered them by first appearance.
BNN – Bayesian neural network
Fermi-LAT – Fermi Large Area Telescope
BLL – BL Lacertae (just one type of blazar)
FSRQ – flat-spectrum radio quasar (just another type of blazar)
DNN – dense neural network
ROC – receiver operating characteristic
AUC – area under curve
References
[1] Jospin et al. “Hands-On Bayesian Neural Networks—A Tutorial for Deep Learning Users.” IEEE Computational Intelligence Magazine, vol. 17, no. 2, 2022, pp. 29-48. IEEE Xplore, 10.1109/MCI.2022.3155327.
[2] Butter et al. “Classification of Fermi-LAT blazars with Bayesian neural networks.” Journal of Cosmology and Astroparticle Physics, 2022. IOPscience, 10.1088/1475-7516/2022/04/023.
[3] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
Week 4: We’ve Officially Entered Orbit
Hello again, and thanks for stopping by for another week of my senior project! Last week, I familiarized myself with the technical aspects of reproducing models and results from prior papers. This week, things get a little more exciting! We’re transitioning away from last week’s experimentation to developing the model for my actual project. Of course, that’s not to say that the work I did before wasn’t related to my project, but it was all as preparation for the big moment.
I know you’re probably dying to hear about this week’s progress, so I won’t keep you waiting any longer with my rambling.
A Small Change in Trajectory
In previous weeks, I was constantly experimenting with smaller datasets. If you’ve been following me along from the very beginning, you’ll remember that I created a smaller dataset with only millisecond and young pulsars. Originally, I thought this would be enough to train my model, but since then, I’ve slightly adjusted my expectations. After extensive discussion with my external advisor, I’ve decided it would be more helpful to distinguish pulsars from other high-energy sources. Rather than simply classifying between millisecond and young pulsars, I’ll be able to contribute more to solving the “candidate selection problem” mentioned in week one’s blog post by actually separating pulsars from other sources.
While this is certainly an important change, it doesn’t significantly affect my project’s trajectory. All of the work we did beforehand still matters. We just have to view things from a different perspective. Instead of using the pulsar dataset we created earlier for training, we can use it as a way to understand the key features of the Fermi Large Area Telescope (Fermi-LAT) catalog[1]. Indeed, working with the smaller pulsar data first gave me a better idea of what features I needed to include in my actual training data.
The Ultimate Dataset
Using my familiarity with the Fermi-LAT catalog, I created the ultimate training dataset this week. If you recall from last week’s blog post, I used a rather small dataset to recreate the results of Butter et al.’s paper on blazar classification[2]. Their dataset included only seven features: seven energy bins with their respective measured energy fluxes. While their model seemed to perform well, this limited dataset might not be enough for my larger classification problem. Since I aim to classify between pulsars and many other types of energetic sources, I’d like a more robust dataset.
To do so, I considered all 21 different types of energetic objects from the Fermi-LAT catalog. Along with pulsars, some of the other sources included are blazars, supernova remnants, galaxies, and molecular clouds. For our purposes, it’s not incredibly important to have detailed knowledge about each source type. It’s enough to know that they’re all high-energy, gamma-ray emitters that contribute to the complex astrophysical background. All of these sources gave us 7,195 data points. This is a pretty substantial jump from the 295 pulsar samples we had originally!
For each of these sources, I gathered a comprehensive list of spectral features that describe each source’s spectral energy distribution (SED). I divided these features into six categories: “identification,” “association,” “flags,” “spectral,” “variability,” and “energy band.” In total, there were 137 feature columns, so I won’t even begin to try listing them all out. Instead, I’ve provided a table with just a few features to give you a general idea:
With all of these features, I was able to construct a large dataset of all sources in the Fermi-LAT catalog. We can see that the dataset includes both the measured values and the uncertainties of the measured values. But of course, it isn’t much use to look directly at the data table itself. After all, it’s a 7,195 x 137 table, so that would be a lot of eye-straining. Instead, I created a couple of counts plots to easily visualize the overall distribution of the data:
After creating these plots, I got a much better idea of how my data was split. For example, looking at panel A, I could see that many of the sources remained unassociated. Meanwhile, most of the associated sources were classified as “BCU,” or active galaxies of uncertain type.
Ready for Training
Parsing through the catalog and mining the data was only the first step–the teeny tiny tip of the iceberg. The larger task was labeling that data for my supervised learning model. This seemed pretty intimidating at first. Deciding how to categorize the data is one of the most crucial steps. It’s what my model is going to be learning from as it trains! To be completely honest, I sat in front of my laptop and stared at my dataset for a couple of minutes before I mustered up the courage to begin the labeling process. However, after discussing with my external advisor, I was able to split up this process up into smaller subtasks.
Rather than defining a single categorical variable to split the data by, I came up with a hierarchy of labels. Splitting the data up this way not only made it easier to label, but it could also help my model learn patterns at different levels of granularity. To complete this step, I used a two-level labeling hierarchy. The higher-level split was between galactic and extragalactic sources. This differentiation separates sources that originate within our own galaxy, the Milky Way, from those that lie far beyond it. It captures fundamentally different physical origins and observational characteristics, which makes it a meaningful first split for the data. Within each of these two categories, I split the sources into their catalogued source type classifications:
Galactic: The galactic group contained msp_like, young_psr_like, and other_galactic objects. As mentioned in earlier posts, millisecond pulsars (MSP) and young pulsars (PSR) are both pulsar objects. They’re primarily galactic, existing throughout our Milky Way galaxy.
Extragalactic: The extragalactic group contained fsrq_like and bll_like sources. Flat-Spectrum Radio Quasars (FSRQ) and BL Lacertae objects (BLL) are both subclasses of blazars. Blazars are extragalactic objects that exist at the center of distant galaxies, not our own.
Along with these labels, I created an extra label for sources that lack confident associations. These included active galaxies of uncertain type, special case objects, and unknown objects. While I did not remove these sources from my dataset, I made sure to flag them as “not trainable” objects. After splitting the data, I obtained 20 extragalactic classes, 23 galactic classes, and 6 not-trainable classes.
As an extra safety check, I engineered an extra column that identified whether sources were identified or had high confidence associations. I used this to separate objects that I considered good enough for training from those that I would only use to evaluate later. In the end, I had 2,877 trainable, high-confidence objects. The other 4,318 were non-trainable objects, whether it was because they had weak associations or were simply unassociated. After excluding all of the non-trainable objects, I created some additional counts plots of the categorized data. They are shown below:
What's Next?
Slowly but surely, we’re developing the pipeline for training our model. Within the next couple of weeks, I’ll be creating feature matrices to be fed into the model and developing the Bayesian neural network. I anticipate a lot of exciting progress, as well as many challenges along the way. I’m more than thrilled to start building and training the model after so many weeks of preparation. Regardless of how it goes, I’ll of course be taking you along with me, through all of the ups and downs. All I can say is buckle up, because things are about to get really interesting!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
Fermi-LAT – Fermi Large Area Telescope
SED – spectral energy distribution
BCU – active galaxy of uncertain type
MSP – millisecond pulsar
PSR – young pulsar
FSRQ – Flat-spectrum Radio Quasars
BLL – BL Lacertae objects
References
[1] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
[2] Butter et al. “Classification of Fermi-LAT blazars with Bayesian neural networks.” Journal of Cosmology and Astroparticle Physics, 2022. IOPscience, 10.1088/1475-7516/2022/04/023.
Week 5: The Future Is in the Features
Hello fellow readers, and welcome back! Isn’t it crazy how we’re almost halfway through the senior project? It felt like it was just last week when we were still getting to know each other! I guess time really does pass a lot quicker when one’s having fun. Talking about fun, I made lots of fun, productive progress this week. I focused on creating the feature matrix that will serve as the main input to my neural network model. So come on! Buckle up and prepare to dive into the feature matrix!
Looking Past the Obvious
Last week, I created a training truth table with all of the sources and whether they were “trainable” or not. If you recall, we defined trainable sources to be those with confident associations. Any sources that weren’t associated or had dubious associations were not considered to be trainable. However, we’re not outright removing them from our data. After all, they’re still valuable datapoints that we’ve simply stored for later.
This training truth table is only the first step. While it categorizes the trainable sources, it’s just a jumbled compilation of sources with no real organization. Using the training table columns as a starting point, I have calculated more useful features. The original training table includes information such as flux measurements across consecutive energy ranges, along with their uncertainties. These merely act as a simple baseline for the spectral energy distribution (SED), but we can derive so many more spectral characteristics and nuances!
Not All Data Is Created Equal
I began this process by splitting my data into a “main energy range” and a “low energy range.” As a reminder, my dataset includes observations of fluxes and energies of sources across a rather large energy range (50 megaelectronvolts (MeV) to 1000 MeV, to be exact). While the data records all of this information for each source, not all of this information is equally reliable.
This is because of something known as a telescope’s point spread function (PSF). A telescope’s PSF describes how the telescope spreads light from a single point source, like a star, into a blurry blob. It essentially characterizes the telescope’s resolution. So why does this matter? I promise I’m not going on a techy tangent for no reason here. This PSF is actually very important to our dataset because it determines which data entries are most sound. When we go below about 300 MeV, the PSF of the Fermi Large Area Telescope (Fermi-LAT) becomes quite broad. It spreads light from the point source over a large area, creating a lower image resolution. This means that the data obtained below the 300 MeV threshold might not be as trustworthy. For this reason, I split the data for each source into low and main energy ranges.
Again, I didn’t completely disregard the low-energy bands. They still contain useful data, but they belong in a nuisance layer rather than the main training layer. We can look at them later to handle systematic noise and background contamination.
Featuring the Features
After creating the energy split, it’s time to move on to some more advanced features! There are two columns in the original training table called Energy_Flux100 and Unc_Energy_Flux100, which describe the energy flux of a source in the 100 MeV to 100 GeV range and its associated uncertainty. I took the logarithm of both these columns. This transformation helps compress the wide range of values in the data, since fluxes can span several orders of magnitude. Working in log space makes these extreme values more manageable, allowing the model to learn patterns more easily. This is especially useful for the uncertainty because it ensures the model interprets errors on a similar scale, rather than being dominated by a few very large values.
I then introduced a derived feature called LP_Index1GeV, which captures the spectral slope of a source specifically at 1 GeV. This value comes from the log-parabola model, which is simply a function commonly used to describe curved gamma-ray spectra. Unlike a simple power law model (which is a straight line in log-log space), a log-parabola allows the slope to change with energy. Put simply, think of a power law as a straight line and a log-parabola as a curve. With this new LP_Index1GeV feature, we can see how the source’s emission behavior changes with energy. It’s much more refined than a single global slope value.
A couple more log calculations later, we finally have the main features for each source. Along with these main, mandatory features, there are a couple of auxiliary features. The log10_LP_EPeak and Frac_Var are included only for sources where these values are finite. They give us information about the peak energy of the source’s emission and the degree of its variability over time, respectively.
Feeling Uncertain?
All of the features I just mentioned represent the essential information in the original Fermi-LAT catalog[1] in a more compact form. Alright, phew! That means we’re done now…right??!?! Actually, not quite! We still have a couple more features to run through before we can safely dock for the night. This next part is arguably one of the most important feature-engineering steps in this entire process. It has to do with the measured uncertainty in the original Fermi-LAT catalog.
What exactly do I mean by this? In the catalog from which we have pulled our data, each energy band does not always come with a clean, confident flux measurement. Instead, the Fermi-LAT pipeline reports either a measured flux or an upper limit, depending on how significant the detection is in that band. If the signal from a source is strong enough relative to the background, it gets recorded as a flux with a reasonable uncertainty. However, that’s not always the case! There are times when the signal may be too weak or indistinguishable from noise. Here, the catalog reports an upper limit. This essentially tells us that the true flux is below this upper limit, although we can’t be sure what that value actually is.
It’s important to note that even in these upper-limit cases, a nominal flux value is still provided, but it comes with a very large uncertainty. This is very helpful to know because it tells us that the distinction between a true flux detection and an upper limit is encoded by how large the uncertainty is relative to the flux itself. We can now leverage this observation to classify each entry.
To put this into action, I’ve constructed a mask for every source in each energy band. The mask categorizes the source’s measurement status into three groups: 2 = reliable flux detection, 1 = upper limit (or nominal flux with a large uncertainty), and 0 = missing entry altogether. For each band, we must first ensure the flux and its uncertainty are finite and positive. It would not make sense for the flux to be negative because it represents the amount of energy received from a source, which can’t be less than zero!
Once we have confirmed that the entries are indeed valid, we can compute the relative error by dividing the uncertainty by the flux. If this value is less than or equal to one, we treat it as a detection; otherwise, we label it as an upper limit. Any entries without a valid positive flux are marked as missing. This mask allows the model to explicitly distinguish between strong and weak measurements. As I did for previous measurements, I’ve taken the log-scaled flux values for both detections and upper limits, while the large uncertainties for upper limits naturally down-weight their influence during training. This allows us to accurately represent the SED without disregarding any important values, however “uncertain.”
Don't Be A Nuisance
For the final step of our feature-engineering process, we have to deal with the low-energy bands. Remember from before that we split the energy bands into main and low energies. This was due to the broad PSF at lower energies, which can make measurements unreliable. Now that we’ve played around with the main energy bands, we can finally deal with these low-energy nuisances.
Because they’re in their own feature block, I can treat them with more care and caution. Just like with the main bands, we should evaluate whether each measurement is a reliable detection or an upper limit by comparing the uncertainty to the flux. From this, we can use the same mask for each band (2 = detection, 1 = upper limit, 0 = missing) and keep both the log-scaled flux and uncertainty values.
Beyond individual flux values, I’ve also engineered a few ratio-based features that compare low-energy flux to mid- and higher-energy fluxes. These ratios help capture how skewed a source’s emission is toward lower energies. It basically gives us a convenient way to describe its spectral shape without relying too heavily on noisy measurements.
There are just a few more supplementary details that provide more useful context for each source. For example, the catalog provides flags that signal potential concerns with a source, such as confusion with nearby sources or issues in the SED-fitting process. These flags are represented by single integers that can be added together to indicate multiple concerns. I’ve unpacked these flag bits into separate binary features, so the model can directly learn from any known data-quality issues. I also included the log of the spectral peak energy where available. You can see the summary of these features in the image below:
Altogether, this approach lets the low-energy data contribute useful information without letting it dominate or distort the main learning process.
Plots, Plots, Plots
Since I've done so much talking with numbers and code blocks, I'd like to send you off with a more friendly, easy-to-understand section. As a last little treat dedicated to you, my esteemed reader, I’ve provided a couple of plots for an easier visualization of some of the mentioned features. Let’s go in chronological order, shall we?
We first looked at the original catalog features, including flux measurements and uncertainties. The plots of Block A below depict kernel density estimates (KDEs) for three key features from the compact catalog: energy flux, uncertainty in energy flux, and average significance. Each of these features is split in three ways: associated vs. unassociated sources, flagged vs. unflagged, and high- vs. low-latitude sources. You can think of these curves as smoothed histograms showing where most sources exist along each feature’s range, with taller peaks meaning more sources crowded around that value.
Some very interesting patterns jump out from here! On the leftmost panel, we can see that unassociated sources (red dashed) consistently skew toward lower energy fluxes and lower significance compared to their associated counterparts. This is as expected because the sources we can't identify tend to be fainter. In the middle column, flagged sources (orange dashed) similarly tend to cluster at lower flux values, hinting that problematic detections are often on the weaker end of the signal spectrum. The high- vs. low-latitude split on the rightmost column is arguably the most visually dramatic. High-latitude sources (purple) form tighter, more defined peaks, while low-latitude sources (cyan dashed) spread out more broadly. This is likely because the low-latitude sources are crowded in the galactic plane, the messiest region. Here, confusion and background noise make detections a lot harder.
Next, we spent a lot of time discussing flux and upper limits in the SED. A few of these log-scaled fluxes and their associated uncertainties are shown across different energy bands in the Block B plots below. We’re again seeing some KDEs that act like continuous histograms, showing the frequency of observed source patterns.
A couple of features are already pretty obvious. The flux distributions (first image) are strikingly narrow and symmetric compared to the plots in Block A, with most sources clustering tightly around similar log-flux values (~-10 to -9). It’s quite surprising to see how similar the associated and unassociated sources look in log-flux space. I actually wasn’t expecting this! But this suggests that at the SED band level, associated and unassociated sources aren't dramatically different in brightness. On the other hand, high- vs. low-latitude splits look a bit more interesting: low-latitude sources (cyan dashed) peak slightly higher in flux and are a bit broader, again reflecting the noisier galactic plane environment.
The uncertainty distributions (second image) look a little more varied. Here, the splits become obviously more pronounced. Unassociated sources also tend to have slightly higher flux uncertainties than associated ones. Flagged sources (orange dashed) and low-latitude sources (cyan dashed) both show much flatter, wider uncertainty curves, meaning their flux measurements are far less precise. This is just as we expect! If we can't pin down what a source is, we probably can't measure it as cleanly either. Together, these two plots highlight the concept of the upper limits we discussed before. While raw flux values may look similar, the reliability of those measurements varies quite a lot.
Hopefully, these plots have given you a better idea of what the underlying data distribution looks like. For one, it’s much easier to understand compared to blocks of code and lines of abstract numbers. These plots are also useful for predicting what patterns our future model might be learning from. Altogether, this week has taught me a lot about looking past the surface-level numbers to really get at the structure of the data. I hope you’ve enjoyed reading as much as I’ve enjoyed preparing this post for you. Until next time, keep looking up!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
SED – spectral energy distribution
PSF – point spread function
Fermi-LAT – Fermi Large Area Telescope
LP – log-parabola
KDE – kernel density estimate
References
[1] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
Week 6: Navigating the Model Galaxy
Hello there, and welcome back! This week is all about exploring baseline machine learning models. Before we jump headfirst into creating our Bayesian neural network, it’s good to see some initial results for comparison. It’s not just mindless, busywork. It’s how we’re going to find out later whether the Bayesian model adds substantial information. To accomplish that goal, I’ve created two strong, baseline models that use the catalog features in a careful, reproducible way. Let’s get better acquainted with them, shall we?
The two machine learning algorithms I’ve chosen for this week are the random forest and CatBoost. Both of these models are built on the premise of decision trees. You can think of decision trees as essentially a series of yes/no questions. It’s quite similar to how we arrive at logical conclusions ourselves! They ask specific questions about the data, narrowing down the possible classifications until a final decision is reached. This makes them a natural fit for teasing apart the subtle differences between types of gamma-ray sources.
This tree-based approach is especially useful for my dataset, where relationships tend to get messy and non-linear. A simpler model would likely not capture all of these complexities. On top of that, these models perform well with minimal preprocessing and handle imbalanced data reasonably well, making them ideal for my project. Don’t worry if you’re still a little confused about the specifics of each algorithm; I’ll get into just that in the following sections. Just keep reading!
A Walk Around the Forest
Ah, yes, the sounds of birds chirping and leaves rustling in this sun-dappled forest. Wait…this isn’t what this post is about, is it? Unfortunately, we’re not talking about actual forests, although the model does branch out in surprisingly similar ways.
As its name suggests, the random forest algorithm is rather like a “forest” of multiple decision trees. More technically, it’s an ensemble of decision trees trained on random subsets of data. At this point, you might start to wonder, “Why go through the effort of training multiple decision trees when you could just do it once?” This is exactly the advantage of using a random forest over a simple decision tree. Instead of relying on a single tree, a random forest averages the results of many trees to make a more stable and reliable prediction. This helps it reduce overfitting and improve generalization to new data, making it generally more accurate and robust[1].
Thankfully, we’re using Python to code, which comes with plenty of machine learning libraries. The scikit-learn library[2] specifically includes many built-in algorithms (called estimators) that we can use to train on our data. It conveniently provides a RandomForestClassifier, which is exactly what we need for this section!
Pre-Training Protocols
Once we’ve loaded all the necessary libraries (including scikit-learn and some other standard utilities), we can define our training sample. We’re using the same training data we curated two weeks ago, and we’re making sure to include only the trainable sources. As a reminder, these are the sources with confident class associations. After filtering out the non-trainable ones, we’re left with 2,877 unique sources that belong to one of five classes: BL Lacertae objects (BLLs), flat-spectrum radio quasars (FSRQs), millisecond pulsars (MSPs), young pulsars (YPs), and another catch-all class of other galactic objects.
Now, we can’t just stick these data points into our classifier and hope for the best. We must first preprocess our data and account for possible sources of bias. The first step is to make our data interpretable to the model. Since the RandomForestClassifier expects the target variable to be an array of integers or floats, we must transform our class labels into numbers. We can use a label encoder, which automatically assigns a number to each class label in alphabetical order. The labels created, along with class counts, are shown in the image below.
At this point, we can see a not-so-tiny problem: there’s a major class imbalance. The biggest class (BLLs) has about 23 times the number of sources as the smallest class (other galactic objects). That’s a big difference! Such a big difference can really affect the training process. If the data were to get split randomly, we could end up with imbalanced class distributions between the training and validation sets. To combat this problem, we can use a strategy called repeated stratified k-fold cross-validation.
Essentially, repeated stratified k-fold cross-validation (wow, that’s a mouthful to say) splits the data into k groups (folds) multiple times, making sure each group keeps the same class balance. This prevents groups from ending up with too many majority-class members or too few minority-class members.
From there, I took things a step further by stratifying not just on the class labels, but also on how bright each source is. I split the data into three “flux terciles:” faint, medium, and bright. I combined that with the 5 class labels to create 15 unique strata. This way, each fold gets a balanced mix of not only source types, but also of faint and bright objects.
Once the data is split, the model trains and tests across different combinations of these folds over several rounds, giving a more reliable sense of performance. In our case, we’re using a 5-fold setup and repeating the whole process 10 times with a fixed random seed. Think of it like reshuffling the same deck of cards over and over while the proportions stay consistent. All of the intermediate predictions are then accumulated and averaged to produce the final result. With that, we get a prediction for the source’s class!
Random Forest Results
To train the data, we’re going to be using the same features we derived from the catalog[3] last week. This includes characteristics like the logarithmic energy fluxes, flux uncertainties, and detection significance, among other things. The RandomForestClassifier digests the data and follows the training process described above. After the model takes a couple of seconds to make predictions, we can investigate its results. Let’s take a detailed look at the results using a numerical table.
From the table above, we can immediately notice a few things. The most striking feature is the poor performance of other galactic objects. However, it’s not surprising, considering that it’s the minority class by a factor of 23. It’s also worth remembering that it’s an umbrella class for all the remaining sources that aren’t categorized into the other classes, which can become difficult to classify. In the coming weeks, I’ll definitely need to dig into better ways to handle this minority class so it doesn’t quietly sabotage the future Bayesian model’s performance.
Other than that one class pulling down the average performance, the other classes seem to be cruising along quite nicely. There’s a solid precision/recall balance, with F1 scores mostly in the 0.74-0.90 range. Of course, the struggling class drags down the macro averages, so we need to be mindful about generalizing the model’s performance. Still, with an impressive ROC-AUC of almost 0.95, it’s clear the model has learned meaningful patterns.
Cats and Boosts
The second model we’re going to discuss today is CatBoost[4]. As I mentioned before, this model is also related to decision trees. Instead of building one big tree like a random forest, CatBoost builds many smaller trees sequentially. Each new tree focuses on correcting the mistakes of the previous ones. It’s a type of gradient boosting algorithm, meaning it gradually improves predictions step by step. What makes CatBoost especially useful is its ability to handle categorical features efficiently and reduce overfitting. This makes it both powerful and stable, a wonderful choice for a baseline model!
To train our CatBoost model, we follow the same preprocessing steps we used for the random forest. We stratify our data for the k-fold cross-validation, and we define our data features. The only difference is that, unlike the RandomForestClassifier, CatBoost does not automatically define class weights for class imbalances. The solution is easy enough; we just need to construct our own class weight dictionary and pass it to the model during training. Just like before, we average the results of ten 5-fold cross-validation processes.
Using these results, we can take a look at the CatBoost model’s overall performance. Once again, we’ll start by deciphering a table of numerical results.
While this table looks almost identical to that of the random forest, I can assure you they’re not the same. Right off the bat, we can see that the model performs poorly on other galactic objects, once again. This is both good news and bad news for us. It tells us that the random forest results were not just a fluke (good!), but it also means that we’ll have to address the minority class so it doesn’t continue to drag down our future models (bad!).
Despite the relatively poor metrics, we still have to give CatBoost credit for performing slightly better than the random forest. Recall from before that the other-galactic-source numbers for the random forest were around 0.17. Meanwhile, the CatBoost average for other galactic objects is around 0.23. It’s small, but it’s an improvement nevertheless. This is likely a result of the model’s internal structure. Since CatBoost builds its model step-by-step, it has a better chance of picking up on subtle patterns in those rare, underrepresented classes that the random forest might otherwise gloss over. This already gives us a better idea of how different models perform on the same data.
Plots and More Plots
Let’s end today’s post with some final plots, because in case you haven’t noticed, I love using visualizations to make sense of things. The first plot is a series of calibration curves for each class. These curves show how well the model’s predicted probabilities match reality. Essentially, if the model is 80% confident, the curves tell us whether it’s actually right about 80% of the time. In a perfectly calibrated model, the curves would lie right along the diagonal line.
From the curves above, we can see some deviations from an ideal scenario. Some classes track the diagonal fairly well, while others tend to drift above or below, indicating underconfidence and overconfidence, respectively. In particular, we can see the minority class exhibiting extreme fluctuation in the bottom left. Overall, the curves suggest that while the models aren’t perfectly calibrated, they’re doing pretty well at turning their predictions into meaningful probabilities. There’s just a little bit of uncertainty spread throughout, which is exactly what we hope to capture with the Bayesian model later.
The next graph, shown below, is a visualization of what we discussed with numerical tables earlier. These two side-by-side heatmaps are confusion matrices for our two models. Confusion matrices represent how often models correctly predict each class (the diagonal) versus how often they confuse one class for another (the off-diagonals). Looking at both plots, we can see the overall structure is very similar. However, we can see once again that they both heavily misclassify the other galactic objects. As we discussed before, CatBoost does show a slight improvement here. It seems to spread predictions a bit more evenly, but the core issue remains. In short, both models are telling similar stories, but still struggling to recognize extra sources in the background.
The final visualization on the menu is a group of precision@k plots. These plots measure the proportion of correctly identified sources within the top k highest-confidence predictions. A higher precision@k indicates that the model’s top predictions are reliably on-target. This is especially relevant to my project because the candidate selection problem depends most heavily on the most confident outputs. Let’s take a look at our very own precision@k plots below.
The three precision@k plots shown above are evaluated at k = 20, 50, and 100. We can see several consistent patterns across the two models. For the well-populated classes, both models achieve pretty high precision across all three thresholds. This performance is largely stable as k grows, which suggests that the high-confidence predictions for these classes remain clean even as the window expands. The millisecond-pulsar class shows a modest but consistent decline as k increases, which might tell us that the models’ lower-ranked confident predictions for this class are slightly noisier. At the very bottom (again) are the other galactic objects, which show virtually no improvement with increasing k. This reinforces the need to study this class separately in the coming weeks.
Looking Ahead
And just like that, we’re done with week six of the senior project! This week was definitely a big step towards developing our Bayesian neural network. I had a lot of fun creating the two baseline models, and I gained some important insights into the data.
In the following weeks, I’ll focus on tackling the minority classes, specifically the other galactic objects and millisecond pulsars. Especially because millisecond pulsars are the core focus of this project, I’d really like to make sure they’re properly accounted for. If you have any suggestions for dealing with such big class imbalances, please leave them in the comments below. I’d love to hear! If not, I’ll see you very soon for some more exciting updates. Farewell…for now!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
BLL – BL Lacertae object
FSRQ – flat-spectrum radio quasar
MSP – millisecond pulsar
YP – young pulsar
References
[1] Kavlakoglu, Eda. “What is random forest?” IBM, 2026, www.ibm.com/think/topics/random-forest. Accessed 1 Apr. 2026.
[2] Pedregosa et al. “Scikit-learn: Machine Learning in Python.” Journal of Machine Learning Research, vol. 12, 2011, pp. 2825–2830. arXiv, doi.org/10.48550/arXiv.1201.0490.
[3] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
[4] Kolli, Aravind. "Understanding CatBoost: The Gradient Boosting Algorithm for Categorical Data." Medium, 13 Feb. 2024, aravindkolli.medium.com/understanding-catboost-the-gradient-boosting-algorithm-for-categorical-data-73ddb200895d. Accessed 3 Apr. 2026.
[5] Zastanceanu, Denis. “Forest path and trees.” Wikimedia Commons, 15 May 2016, commons.wikimedia.org/wiki/File:Forest_path_and_trees.jpg. Public Domain.
Week 7: Into the Bayesian-verse
Hello, hello, and welcome back! We’re seven weeks in, and I’m very excited to announce that we’ve finally begun building our Bayesian neural network. This is the moment we’ve all been waiting for, isn’t it? I won’t keep you waiting any longer, since I can literally feel your excitement pulsing through the screen. Without further ado, let’s get straight into this Bayesian business!
Let's Rewind
Just as a small reminder, we trained two baseline models last week: random forest and CatBoost. While both models performed well overall, we noticed they had trouble with the minority class–the other galactic objects. The evaluation metrics for that class were alarmingly low compared to the other classes. However, after discussing it with my advisor, I’ve concluded that it’s not a major concern. At the end of the day, the other-galactic class is not a coherent physical class, and both models correctly reflected that uncertainty. Since my project doesn’t rely on other galactic objects, it seems reasonable to leave this aside for now. All is well as long as my models can accurately differentiate between the other specific classes (which they seem to be doing quite well at the moment).
With that being said, I’ve decided to move on to the core Bayesian part of my project. Of course, our results from last week will serve as fundamental comparisons. We’ll specifically see how the Bayesian models capture model uncertainty better (or worse!) than our baseline models.
Starting Simple
This week, I constructed a simplified version of a Bayesian neural network (BNN). It’s not exactly a proper BNN, but it closely resembles one. Let me explain. What I create this week is essentially a Bayesian multilayer perceptron (MLP). The only difference is that, instead of giving each network weight a probability distribution (which is what a “true” BNN would do), I follow a much simpler idea. This idea comes from the paper “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning” by Yarin Gal and Zoubin Ghahramani[1]. The paper introduces the idea that we can use dropout as a form of approximate Bayesian inference. As the term suggests, “dropout” is when we randomly drop neurons during training. By doing this, we’re implicitly sampling from a distribution over many possible neural networks.
I’m choosing to take this route first because a Bayesian MLP is much simpler to implement than a traditional, full-scale BNN. Developing a BNN is no small task, and it involves very advanced machinery like variational inference and reparameterization tricks. So instead of diving straight into this complexity, I can start with a familiar neural network and add in uncertainty. This makes the transition from standard models much more intuitive.
Just as importantly, this approach sets up a useful comparison for next week. By starting with a Bayesian MLP, I’ll have a lightweight approximation for how much uncertainty estimation can improve performance and interpretability. Then, when I move to a proper BNN with learned weight distributions, I can directly evaluate whether that extra complexity brings any additional benefits. And you’ll find later that even at this simplified stage, things already get more complicated compared to the baseline models of last week. Introducing dropout and Bayesian inference adds a whole new layer of depth. This will all show up clearly in the results I’ll dive into later in this post. For now, let’s move on to the model architecture.
An Architecture Lesson
The architecture starts off looking pretty familiar. It’s just a standard neural network with fully connected layers, ReLU activations, batch normalization, and dropout after each hidden layer. However, that dropout is doing much more than usual. Instead of just being a strategy to prevent overfitting, it’s the core step that turns this into a Bayesian model. Every time data flows through the network, different neurons are randomly “dropped,” or ignored. This means the model is never the same across different passes.
Training the model follows a standard pipeline, but with a few important twists. We first split the data using a cross-validation scheme, where each source appears in multiple validation folds across different repeats. Since we’re using the same repeated stratified k-fold cross-validation technique as we did last week, we can reuse the same fold splits for consistency. Thankfully, I’ve saved all of the fold assignments, which makes this task much simpler.
For any given fold, the model trains on the training split and evaluates on the validation split using standard cross-entropy loss. During this phase, dropout behaves normally and helps the model generalize better. However, there’s an extra step in the middle, known as early stopping. This is where we temporarily turn dropout off and switch the model to evaluation mode to compute a stable validation loss. If that validation loss stops minimizing, we can stop training early to prevent overfitting. This is an especially popular approach when training MLPs to control model complexity[2].
To make the actual predictions, we use the familiar repeated stratified k-fold cross-validation scheme. This setup also enables out-of-fold (OOF) predictions, which are key for fair model evaluation. For each fold, the model predicts only on the validation data it hasn’t seen during training. These predictions are then accumulated across all folds to form a complete set of OOF predictions for the entire dataset. It mimics how the model would behave on new, unseen data. In our case, it sets up a direct comparison with the baseline models from last week, since they use the exact same fold structure.
Finally, the real magic happens at inference time with Monte Carlo (MC) dropout. This is the concept introduced by Gal and Ghahramani, which I mentioned previously. Instead of turning dropout off during inference, the model is kept in training mode, with dropout staying active. For each input, the model performs many forward passes (we’re using 50 today), each time with a slightly different subset of neurons active. This produces a collection of predictions, which are then averaged to get the final probability. The spread of those predictions quantifies the uncertainty we’ve been looking for for so long.
Hyperparameter Tuning
Another step we have to look at before training is hyperparameter tuning. This is definitely the least exciting part, but it’s also one of the most important aspects of the model’s architecture. The difference between two sets of hyperparameters could really make or break our model. To accomplish our goal, we’re using an open-source hyperparameter framework called Optuna. Optuna actually uses Bayesian optimization itself to search the hyperparameter space. Rather than exhaustively trying every combination like grid search, it builds a probabilistic model of which regions of the search space are likely to yield good results, focusing subsequent trials there.
For my Bayesian MLP, I run 50 trials searched over five hyperparameters: network depth, layer width, dropout rate, learning rate, and batch size. Each trial trains on the MLP for up to 60 epochs with cosine annealing, which is a technique that gradually reduces the learning rate to help the model reach minima in late training.
The left panel of the plot below shows the optimization history, which tracks the validation negative-log-likelihood (NLL) across trials. NLL measures how confidently and correctly the model assigns probability to the true class labels, and we want to minimize this value. Clearly, we can see the best-so-far NLL drops sharply in the first ~10 trials, ultimately bottoming out at NLL = 0.0750. Meanwhile, the right panel shows a parameter sensitivity chart. We can see that batch size and learning rate are by far the most influential parameters.
After complete tuning, the optimal configuration reflects those sensitivities quite well. We end up with a relatively wide, two-layer network (2 x 512 units) with aggressive dropout to regularize it. We also get a moderate learning rate of ~0.0024 and a small batch size of 64. The high dropout rate paired with a wide network architecture is commonly seen, as the width provides expressivity while dropout prevents co-adaptation between neurons. The small batch size, which is the most sensitive in our case, likely helps by introducing gradient noise that aids generalization. All of these optimal parameters are summarized in the table below.
Parameter | Optimal Value
# of layers | 2
Layer width | 512
Dropout rate | 0.47
Learning rate | 0.002
Batch size | 64
Creating a Hierarchy
Now that we’ve got all the prep work done, we can finally get to training our model. Rather than throwing all five classes at a single neural network and hoping for the best, we’re going to take a divide-and-conquer approach. Instead of giving a single neural network a massive classification task, we’re going to create a hierarchical classifier. You can think of it essentially as a two-stage decision tree made of neural networks. The first network, which we’ll call Level 1, makes the first split by asking the big-picture question: is this source extragalactic or galactic? Once that broad decision is made, two separate Level 2 networks each handle the finer-grained classification within their respective branch. One sub-classifier distinguishes between the two extragalactic types (BLL-like vs. FSRQ-like), and the other separates the three galactic types (MSP-like, other galactic, and YP-like).
Then, we can compute the final probability for a source as the product of the two predicted probabilities. Let’s say we’re finding the class probability for a young pulsar. We’d take the final probably as P(galactic) × P(young pulsar | galactic). This hierarchical structure simplifies each network’s task, which tends to produce better-calibrated probabilities.
One detail worth highlighting is how uncertainty propagates through the hierarchy. Because each branch produces not just a mean probability but also a standard deviation (courtesy of the MC dropout that makes these networks “Bayesian”), we have to compute the final combined uncertainty for each class. To do this, we use the product rule for error propagation, which looks like
σcombined ≈ σL1 × PL2 + PL1 × σL2.
In plain English, the overall uncertainty in a final class prediction reflects both how unsure the Level 1 network was and how unsure the Level 2 network was about the fine class within that branch. That is ultimately how we capture the uncertainty using this MC-dropout Bayesian MLP model.
Temperature Scaling
As a final step before we arrive at our results, we need to perform something known as temperature scaling. What does model-training have to do with temperature? Everything, actually. The name comes from the Boltzmann-Gibbs distribution in statistical physics[3], where temperature controls how spread out probabilities are across energy states. So in a physics-meets-machine-learning way, temperature scaling makes a model’s predicted probabilities more trustworthy. It specifically calibrates the model to address over- or underconfidence without changing the underlying model itself.
The process is actually quite simple and intuitive. In a normal neural network, the last layer produces raw scores called logits. It then applies a softmax function to turn these logits into probabilities that sum up to one. In temperature scaling, we insert an additional scalar T to divide all the logits before applying the softmax function.
Pi = eyi / TΣk=1n eyk / T
Following some basic math principles, we can see that a larger T (>1) squashes the logits closer together. This creates a flatter distribution and takes mass away from the dominant classes, making it a good fix for overconfidence. On the other hand, a smaller T (<1) amplifies the difference between logits, lending even more confidence to its highest class. This makes it a good fix for underconfidence.
In my case, after training the Bayesian MLP, I apply temperature scaling on a validation set to fine-tune how well those predicted probabilities align with reality. The model finds T by trying out different values on a held-out set with known true labels to see which value of T creates probabilities that are most consistent with the actual outcomes. We fit T on all 50 OOF held-out set predictions and average across them to get our final T. The graph below shows the distributions of T values for different folds, along with the optimal, fold-averaged value we end up with.
Finally, the Results
After training on each fold, we can accumulate and average the predictions to get our final results. Since we have over 2,800 sources, it wouldn’t be very useful to inspect individual predictions. Rather, let’s explore some evaluation metrics. We’re specifically going to look at the accuracy, F1 score, area under the receiver operating characteristic (ROC) curve, NLL, Brier score, and expected calibration error (ECE). Don’t worry if these sound complex and unfamiliar; I’ll explain each one below!
Accuracy
Accuracy is probably the most basic evaluation metric, and it should sound pretty familiar to most. It’s calculated the same way as any accuracy score–dividing the number of correct predictions by the total number of predictions made. However, because it’s the simplest score, it’s also not as comprehensive or trustworthy. That’s why we have to consider a few other metrics to get a better picture.
F1
Unfortunately for my Formula 1 fans, the F1 score is not related to racing. In machine learning, the F1 score summarizes how well our classifier correctly identifies each class. It balances two competing concerns: precision (of all the sources the model said belonged to a class, how many actually did?) and recall (of all the sources that truly belonged to a class, how many did the model successfully identify?) The F1 score takes their harmonic mean: F1 = 2 × (Precision × Recall) / (Precision + Recall). A score of 1.0 means perfect precision and recall, while a score of 0 means the model is doing about as well as a coin flip. In most cases, the higher the better.
ROC-AUC
This metric comes in two parts: the ROC curve and the area under the curve (AUC). We often refer to it simply as ROC-AUC. The first component is the ROC curve, which is the actual plot that depicts how well a binary classifier model distinguishes between two classes across all possible classification thresholds. A well-calibrated model would have a curve that tightly hugs the top-left corner, while a random model would have a curve that lies on the diagonal. Meanwhile, the AUC captures, as the name suggests, the area under this curve. It ranges from 0 to 1, with higher scores indicating better performance.
Negative-Log-Likelihood
I briefly mentioned this metric when we discussed hyperparameter tuning, but let’s look at it in a little more detail. NLL is basically a loss function measuring how well a model’s predicted probability distribution matches true labels. If, for example, a model assigns a 90% probability to the correct class for one source and only 10% for another, NLL will penalize that second prediction far more harshly than the first. Mathematically, it does this by taking the logarithm of the predicted probability for the true class. This turns high-confidence correct predictions into small penalties, and low-confidence or wrong predictions into large ones. It then negates it so that lower is better. Unlike our previous metrics, we want the NLL to be lower.
Brier Score
Of all the probabilistic metrics, the Brier score might be the most intuitive. It’s simply the mean squared error between the model’s predicted probabilities for each class and the actual outcome. Like NLL, the Brier score rewards calibrated confidence rather than just correct predictions, but it’s a bit gentler since it uses squared errors rather than logarithms. Again, we want a lower Brier score.
Expected Calibration Error
The final metric on the list is the ECE. ECE is calculated by first grouping predictions into confidence bins and measuring the gap between average confidence and actual accuracy within each bin. The ECE is then the weighted average of those gaps across all bins. Like the two previous metrics, a low ECE is generally more desirable, and it indicates that our model is well-calibrated.
CatBoost vs. BayesianMLP
Using the metrics we just went over, we can compare the performance of last week’s CatBoost model to that of our Bayesian MLP. I’ve created a quick table below to summarize the evaluation metrics for both models:
Accuracy | F1 | ROC-AUC | NLL | Brier | ECE
CatBoost | 0.8394 | 0.8413 | 0.9455 | 0.4527 | 0.0460 | 0.0299
BMLP (temp-scaled) | 0.8300 | 0.8188 | 0.9400 | 0.4844 | 0.0499 | 0.0118
Looking at the table, we’ll notice that CatBoost, on average, seems to perform better than the Bayesian MLP. This suggests that, as of right now, CatBoost is the stronger discriminative classifier overall. The one surprising exception is the ECE: the temperature-scaled Bayesian MLP clearly outperforms CatBoost here. This might mean that while the MLP may be slightly less accurate in raw terms, its probability outputs are better calibrated. I’ll have to see whether that calibration advantage is worth the hit in raw accuracy. For our particular use case, knowing how confident the model is matters as much as the prediction itself, so the MLP’s superior ECE is a promising sign.
Okay, I think we’ve had enough with numbers for today. Let’s move away from numbers to look at a plot comparing the calibration curves of the two models. The curves below give us a more nuanced picture.
We can see that, for the two most well-represented classes (bll_like and fsrq_like), the Bayesian MLP follows the perfect diagonal much more closely. On the other hand, the Bayesian MLP is more erratic in the other, smaller classes. It’s interesting to note that these other three classes are all from the galactic branch, which is considerably more crowded than the extragalactic branch. That inherent data noise could partially explain the Bayesian MLP’s difficulty in identifying those classes. The other_galactic class is essentially a disaster for both models, which is pretty expected given how few training examples exist. I’ll continue looking into this problem and discussing it with my advisors to decide whether and how to deal with this class.
Looking Ahead
Now that we've gotten familiar with the Bayesian architecture through a simplified MLP, I think it's time we tackle the big challenge. Next week, I'll begin developing the proper BNN with all of its probabilistic layers and complexities. I'll also be considering the problems we faced this week, including the Bayesian MLP's lower performance and the other_galactic minority class.
To address the MLP's current weaknesses, I'll likely reconsider its shallow architecture and/or tune its hyperparameters. Even though it performed slightly worse than CatBoost, we need to do a bit more tuning before we can conclude that the Bayesian approach is wrong for this problem.
Overall, we made some considerable progress this week! I will check back next week to provide more exciting updates about the BNN. Goodbye for now, my fellow reader!
Acronyms
Here’s a list of all the acronyms used throughout this paper (there are quite a few this week). I ordered them in terms of first appearance:
BNN – Bayesian neural network
MLP – multilayer perceptron
OOF – out-of-fold
MC – Monte Carlo
NLL – negative-log-likelihood
BLL – BL Lacertae object
FSRQ – flat-spectrum radio quasar
MSP – millisecond pulsar
YP – young pulsar
ROC – receiver operating characteristic
ECE – expected calibration error
References
[1] Gal, Yarin, and Zoubin Ghahramani. “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.” Journal of Machine Learning Research, vol. 48, 2016, pp. 1050-1059. doi.org/10.48550/arXiv.1506.02142.
[2] Vehtari, Aki, and Jouko Lampinen. “Bayesian MLP neural networks for image analysis.” Pattern Recognition Letters, vol. 21, no. 13-14, 2000, pp. 1183-1191. ScienceDirect, doi.org/10.1016/S0167-8655(00)00080-5.
[3] Mattei, Pierre-Alexandre, and Bruno Loureiro. “The Well-Tempered Classifier: Some Elementary Properties of Temperature Scaling.” 2026. arXiv, doi.org/10.48550/arXiv.2602.14862.
Week 8: Houston, We Have Our Model
Hello there, and thanks for stopping by my senior project again! We’ve officially reached week eight of this long journey, and this week is all about the Bayesian neural network (BNN) we’ve both been dying to hear about. Yes, you read that right. After so much waiting and preparation, we’ve finally gotten to the real meat of this project. But wait…didn’t we already create a BNN last week? Well, the answer is both yes and no.
In case you haven’t read last week’s post (stop and read it!) or have simply forgotten some of the content, I’ll give you a quick recap. Last week, we developed a kind of simplified BNN. We specifically created a Bayesian multilayer perceptron (MLP), using dropout at inference time to simulate different weight distributions. Since I’ve already described so many of the technical concepts in last week’s post, I won’t bother you with them again here. All you need to know is that it was a stripped-down version of a traditional BNN that we’ll use for comparison this week. If you’d like to read more about the Bayesian MLP we used, you can read last week’s post or Gal and Ghahramani’s paper[1], which inspired our Bayesian approximation.
Distributions, Not Points
Enough with my babbling! Let’s jump straight into the BNN, which is what I’m sure you are actually here for. Before we discuss our training results, it’s a good idea to familiarize ourselves with the model architecture. I won’t be talking about the math or statistics behind the model today, but if you’d like to brush up on those concepts, then I encourage you to read this previous post: Week 3: A Bayesian Test Flight.
Last week, we used traditional neural network (NN) layers with point values for weights. This week, we’re using the pre-built DenseVariational layer from TensorFlow Probability (TFP). TFP is a very useful, open-source Python library built on TensorFlow for probabilistic designs, which is exactly what our model needs! Essentially, DenseVariational layers are Bayesian dense layers that use variational inference. They allow the model to learn a probability distribution (posterior) over the weights and biases rather than point estimates. We can then sample from these distributions to capture the model uncertainty by effectively forming an ensemble of different networks.
The architecture itself is pretty straightforward, but it’s powerful nonetheless. We first normalize our inputs, then we pass them through a series of hidden layers, each applying a nonlinear activation like ReLU. We use classic cross-entropy for the loss function, but we also add another term: the Kullback-Leibler (KL) divergence. The KL divergence is what actually makes this model “Bayesian.”
During training, not only is the model trying to fit the data, but it’s also trying to stay close to our original assumptions, which we call the prior. The KL divergence measures how different the learned weight distributions (the posterior) are from the prior by adding a penalty if they stray too much. We scale this penalty by 1/N, where N is the number of training samples. This ensures that it doesn’t overpower the data fit. Combining accuracy with original beliefs, we create a more cautious, mindful model that maintains a good level of uncertainty.
The Training Arc
Now that we’ve got a general idea of the model architecture, we can start discussing our training loop. If you’ve been following along over the past few weeks, this part will look quite familiar. That’s because I’ve been keeping the data preprocessing and training loop pretty similar across the different models to allow for consistent comparisons.
As always, we start by loading and preparing our data for the model. We’re using the same core features from the feature matrix we created several weeks ago to train on our 2,877 sources. Just as we’ve been doing for the past few baseline models, we’re going to use repeated stratified k-fold cross-validation for training. To do this, we first map each source object to its original row index in the training set and use the stored global fold labels to identify which samples belong to each validation split. This reconstruction process preserves the same 50-fold (10 repeats x 5 folds per repeat) structure used previously.
Moving on to the next step, we tune the hyperparameters using the same Optuna framework as last week. This time, we’re no longer using dropout at inference, so we no longer need the dropout rate as a hyperparameter. That leaves us with the number of layers, layer width, batch size, and learning rate to optimize. After tuning the hyperparameters, we get these optimized values, summarized in the table and graph below:
Parameter | Optimal Value
# of layers | 1
Layer width | 256
Learning rate | 0.0097
Batch size | 256
We use the optimal hyperparameter values to develop our model, which follows a training loop similar to last week. For each fold, we use a hierarchical BNN that first predicts whether a source is extragalactic or galactic, and then distinguishes between the specific subclasses. During validation, each source passes through all three models.
At the end, we combine the predictions in a probabilistic way, with the top-level prediction acting as a weighting factor that determines how much influence each branch should have. We repeat this process across all folds and average the predictions so that each source benefits from multiple independent evaluations. After we’ve accumulated all the predictions, we apply temperature scaling to better calibrate the model’s over- or underconfidence. That, more or less, is the complete training process, compressed into a couple of paragraphs.
Model Showdown
It’s time for the model showdown! Now that we’ve got our BNN up and running, we can finally compare it to the baseline results from previous weeks. Let’s specifically look at how the BNN compares to the random forest, CatBoost, and Bayesian MLP models. I’ve created a quick table to summarize all of the evaluation metrics below.
Accuracy | F1 | ROC-AUC | NLL | Brier | ECE
RandomForest | 0.8498 | 0.8416 | 0.9475 | 0.4472 | 0.0439 | 0.0352
CatBoost | 0.8394 | 0.8413 | 0.9455 | 0.4527 | 0.0460 | 0.0299
BMLP (temp-scaled) | 0.8300 | 0.8188 | 0.9400 | 0.4844 | 0.0499 | 0.0118
BNN (temp-scaled) | 0.7588 | 0.7491 | 0.9193 | 0.6061 | 0.0662 | 0.0388
If some of these metrics sound unfamiliar, you can go back to last week’s post to read all about them! The main thing to remember is that we generally want accuracy, F1, and ROC-AUC to be closer to 1, but we want NLL, Brier, and ECE to be closer to 0. Just by looking at the table, we can see that neither Bayesian model performs as well as random forest or CatBoost.
Although the results seem discouraging, they are pretty much what I expected. BNNs tend to be inherently more cautious. They can trade some of the raw accuracy for better uncertainty estimation, which can hurt point predictions on finite datasets. On top of that, they’re definitely harder to train and tune, especially with limited data. Meanwhile, tree-based models like random forest and CatBoost are extremely efficient at capturing structured patterns with less sensitivity to hyperparameters, giving them a practical advantage here. However, don’t get disheartened, my brave reader. This is simply a benchmark result telling us to spend more time reinventing and refining our BNN.
Battle of the Bayesians
Now, if we take a closer look at the table, we can see that the Bayesian MLP seems to outperform the BNN on all fronts. This is likely due to the Bayesian MLP’s simplicity. After all, it’s not a traditional BNN in the sense that it uses dropout to act as a lightweight Bayesian approximation. It’s much less complicated than the full BNN, which must incorporate KL-divergence to balance data fit with prior beliefs. This adds more complexity and noise during training, making convergence harder.
Let’s observe each Bayesian model’s uncertainty plots to get a better idea of where the BNN is struggling. The plots below show the models’ uncertainties against their prediction confidence. In a well-calibrated model, we’d expect the uncertainty to decrease with higher confidence, or a negative slope. Based on that knowledge, we can clearly see that the BNN struggles to tame the uncertainty in confident predictions. While the Bayesian MLP (left panel) shows a clean downward trend, the BNN’s (right panel) uncertainty remains relatively high and scattered, even at high confidence levels. This suggests that the BNN is less well-calibrated, which is what we have to focus on in the following weeks.
Looking Ahead
Now that we’ve seen how our models compare with each other, I have a better idea of how we need to improve our BNN. One of the main issues is that the BNN is dominated by variational inference and KL-divergence, which can overpower the patterns in the data. We have to balance the mean-field Gaussian prior and relatively strong KL weight term so that the model doesn’t constantly get pulled back toward the prior belief.
All of this comes to show how complex BNNs are. They’re very sensitive to approximation choices and tuning, so the slightest mistake can break down the uncertainty drastically. In the following weeks, I’ll be working to fine-tune our BNN so that it’s able to balance knowledge from both the prior and the posterior. It will probably call for a lot of conservative measures, so we’ll experiment together over the next few weeks. If you’re familiar with developing BNNs, I’d love to hear some suggestions in the comments. Otherwise, I’ll see you again in my next post!
Acronyms
As usual, here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
(B)NN – (Bayesian) neural network
MLP – multilayer perceptron
TFP – TensorFlow Probability
KL – Kullback-Leibler
References
[1] Gal, Yarin, and Zoubin Ghahramani. “Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning.” Journal of Machine Learning Research, vol. 48, 2016, pp. 1050-1059. doi.org/10.48550/arXiv.1506.02142.
Week 9: Hope on the Horizon
Hello, hello! Thanks for stopping by again! This week is special in more than one way. For one, I’m typing this out at an airport gate as I await my flight to China. The plane that will soar high into the sky and bring me there is a perfect metaphor for the progress I made this week. If you’ve been following along, you’ll know that the Bayesian models have experienced a few hiccups. But this week, things seem to be taking a turn for the better!
We Can Do Better
After discussing with my external advisor, I was able to pinpoint a couple of areas for improvement in my Bayesian neural network (BNN). The main problems I noticed were related to the input features and training loop choices. I’ll break down each tweak individually in the following subsections, so keep reading to find out!
Never Enough Data
If you recall, I’d created two subsections in my feature matrix. The first one, which we’ll call Block A, contains seven compact features from the Fermi Large Area Telescope (Fermi-LAT) catalog[1]. These features describe the overall properties of each source, such as its spectral shape and variability. The second one, which we’ll call Block B, consists of spectral energy distribution (SED) data. It shows the specific logarithmic flux values of each source across high-energy ranges.
Originally, I included only Block-A features for simplicity because they were the most comprehensive core features. This week, I’ve decided to include the Block-B features as well for a more detailed representation of the data. Here’s a list of both Block-A and Block-B features we’re feeding into our network this week:
Block-A Features | Block-B Features
log10_Energy_Flux100 | flux_300_1000
log10_Unc_Energy_Flux100 | unc_300_1000
log10_Signif_Avg | flux_1_3
LP_index1GeV | unc_1_3
LP_beta | flux_3_10
LP_SigCurv | unc_3_10
log10_Variability_Index | flux_10_100
unc_10_100
flux_100_1000
unc_100_1000
As you can see, the Block-B features are much more extensive. They provide a much better picture of the SED shape across energy ranges, but I didn’t include them initially because there are quite a few missing entries for these features. Nevertheless, they encode important information that could benefit our model. By using these features only when they are available, we can hopefully gain some extra points for our model.
Keep It Simple
The next problem is related to the hyperparameters. Last week, I used Optuna to tune my hyperparameters. I basically used one of the cross-validation folds to tune the hyperparameters and applied those optimal values across all folds. While it seems like an efficient way to obtain optimal hyperparameters, it’s not very stable. Especially because I’m tackling three separate subproblems in my hierarchical BNN, I think it’d be better to stick with conservative default parameters. Once the repaired model is stable, I can consider adding back nested tuning. For now, let’s keep things simple.
A Little Remodeling
The final changes concern the model architecture itself. I’m making some small but important changes to how I structure the input data, compute the uncertainty, and fit the temperature scaling.
Data Normalization
In my original code, I added a BatchNormalization to the BNN input path. BatchNormalization is a popular data preprocessing step that normalizes input features across the batch dimension. It splits the data into multiple batches of samples. It then computes the mean and variance of each feature across all samples in that batch to normalize. The problem with this is that, given my relatively smaller dataset, the batches are small. With small tabular batches, training and inference behavior can become noisier.
Instead of normalizing across a batch of samples, we’re using per-sample LayerNormalization. Here, the normalization happens across the feature dimension for each sample independently. That means no more batch dependency, and no more batch dependency means that the behavior is the same at training and inference. This ensures that the normalization layer behaves consistently, making the entire process much more consistent.
Uncertainty
The next step is to fix how we calculate the predictive uncertainty. Last week, I made the mistake of combining the standard deviations for each of my three submodels in quadrature. I assumed that these standard deviations were independent and that they would propagate linearly into the final output. However, both of these assumptions are wrong. First of all, the uncertainties for the submodels are clearly related to each other because they depend on the same input. Second of all, the propagation is nonlinear. Let’s say we want to calculate the final probability for bll_like. We would get
P(bll_like)=P(extragalactic) x P(bll_like | extragalactic)
This is clearly a product of two random variables. Because we have a product of uncertain quantities, the final uncertainty is not simply the quadrature sum of the input uncertainties. The actual uncertainty in the product depends on both values and their joint distribution in a nonlinear way. It’s certainly more complex, but it’s a change we have to make if we want our model to improve.
To solve this problem, we’re computing the joint uncertainty through all three submodels. Rather than running each model’s passes separately and then combining the standard deviations, we simultaneously sample from the weight posteriors of all three at once. We get a complete set of weights that we can run through the full model hierarchy. This gives us a completely 5-class probability vector rather than a derived approximation.
Temperature Scaling
Last but not least, we’re changing how we perform temperature scaling. Last week, I fit the temperature on all out-of-fold (OOF) predictions. While that sounds reasonable, it creates a little bit of leakage. When we fit temperature on the pooled predictions, some of those very same predictions are the ones we’re later calling our “evaluation results.” This might make some of the evaluation metrics better than they’d actually be on new, unseen sources.
For a more accurate temperature fit, we’re giving temperature scaling its own dataset inside every fold. Essentially, we create a three-way split, where one chunk of data trains the model, one chunk finds the best temperature, and one chunk tests. Each fold also gets its own temperature instead of one global value shared across everything. The resulting metrics will hopefully be more reliable, as they’re genuinely measuring performance on data that was hidden from every part of the pipeline.
Did It Work?
So we’ve talked about all the fixes, but the real question is: do they actually improve our model? Thankfully, it seems like they do! The new BNN, with all of the aforementioned changes, has improved on many fronts. Just looking at the comparison table below, we can see that the repaired BNN has better values for all evaluation metrics. The accuracy, F1, and ROC-AUC are all higher, while the NLL, Brier score, and ECE are all lower.
Accuracy | F1 | ROC-AUC | NLL | Brier | ECE
CatBoost | 0.8394 | 0.8413 | 0.9455 | 0.4527 | 0.0460 | 0.0299
BNN | 0.7588 | 0.7491 | 0.9193 | 0.6061 | 0.0662 | 0.0388
Repaired BNN | 0.7810 | 0.7762 | 0.9328 | 0.5450 | 0.0570 | 0.0285
Although these improvements are small, they’re an optimistic step toward the right direction. We now have a better idea of what the model architecture might need more or less of. If we look at this next plot, we can see that the fix to the uncertainty calculation is a really significant choice.
Immediately, we notice that the repaired BNN’s uncertainty (right panel) has a much more obvious negative correlation with confidence. In other words, the uncertainty of the repaired BNN is decreasing as prediction probabilities get higher, which is exactly what we want. The plots’ axes reveal another small detail. If you look at the y-axis of each plot, you’ll see that the maximum uncertainty value for the initial BNN (left panel) is 0.66, while that of the repaired BNN is only 0.14. So not only are the repaired BNN’s uncertainties less scattered, but they’re also smaller overall! That’s great news because it means our new computation method more accurately reflects the model’s uncertainty.
Looking Ahead
We certainly made some very meaningful progress this week, but it’s not all sunshines and rainbows yet. Looking back at the comparative uncertainty plots, you’ll notice something weird going on around the 0.3-0.4 probabilities for the repaired BNN. There seems to be an area of exceptionally small uncertainties, even for wrong predictions at low confidence. This is one of the missing pieces I’ll have to look into next week.
Along with this small caveat, I’m also planning to investigate the effects of additional small changes. For example, I might consider re-examining the imputation method for Block-B features or seeing whether nested hyperparameter tuning can be safely reintroduced. Since the current model relies on conservative defaults, there is likely still some performance left that could be revealed by a more robust tuning process. By refining how we handle these details, I hope to close the gap even further between our BNN and the CatBoost baseline. All in all, it's been quite a productive week, and I can't wait to see what the results look like once I'm back on solid ground. I'll catch you all in the next update!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
BNN – Bayesian neural network
Fermi-LAT – Fermi Large Area Telescope
SED – spectral energy distribution
OOF – out-of-fold
References
[1] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
Week 10: Rewriting the Stars and My Model
Progress is never linear, and this week, I’ve really started to understand what that means. With only a couple of weeks left, you might want to assume that everything is slowing down and coming to an end. However, that’s not quite the case. We’re still going full steam ahead. This week, I’ve worked with my external advisor to devise new changes to improve the model’s evaluation metrics. Make sure to buckle up, because we’re in for a bumpy ride as we rewrite the model!
Room for Improvement
Last week, we saw small improvements in our Bayesian neural network (BNN) after making a small fix to the algorithm. While the BNN is heading in the right direction, it still performs worse than the uncalibrated baseline models. Ideally, we want to see the BNN performing at least as well as the two baselines, if not better. To this extent, there’s still some room for improvement. In the following sections, I’ll break down the major areas for improvement I’ll be working on over the coming weeks.
Break to Build
You’ve got to sacrifice something to gain something, right? Perhaps the biggest pivot this week is how we’re reconstructing the entire model architecture. Previously, we’d been working with a hierarchical model, where we would train three separate BNNs that each worked on a small classification task. The top-level BNN would identify galactic and extragalactic classes, while the two lower-level BNNs would distinguish specific object types within each class. While it seems like a logical approach, it’s extremely computationally expensive. Just think about how much extra time and space are needed to train three separate BNNs!
To subvert this problem, I’ve switched to using a single, shared-trunk hierarchical model. Let’s break this down a little. By “single,” I mean that we will only be training one BNN. If all goes to plan, that would save us from a lot of computational complexity. Moving on to the “shared-trunk” part, we’re essentially using one big, common trunk for the neural network that acts as the main feature extractor. It takes in all the data and input features to produce a shared internal representation of each source. Once the trunk’s gotten a general idea of the data, it splits into three separate branches. This is where the “hierarchical” part comes in. The three branches serve the same functions as the original three BNNs, but instead of being separate models, they’re now integrated as the final layers of a single network.
In this new model, only the final three branches are Bayesian. The main neural network body that does most of the feature extraction is a regular deterministic model with point-value weights. Prior research in tabular deep learning has shown that simple residual-network-like (ResNet) backbones are strong baselines for tabular problems[1]. ResNets handle deep feature interactions well due to their skip connections, which is exactly what we need for our classification problem. Following these findings, we use a ResNet-like structure for the deterministic trunk. Only then do we attach the three Bayesian branches for fine-grain identification tasks. This further simplifies the model by introducing the Bayesian complexities only in the final layers.
No Data Is Data
Now that we’ve got the biggest change out of the way, we can move on to finer improvements. One problem I noticed is the missing entries in the second block of the feature matrix. This, which we called Block B last week, represents the high-energy spectral energy distribution (SED) measurements of each source. Last week, we simply used median imputation to fill in the missing entries, but didn’t actually expose the missing data. By doing this, we were losing valuable information because the fact that a value is missing carries its own significance.
In our new network, we create an additional layer of information to record the absence of data. We still keep the imputed median values so that our neural network learns properly, but we also add binary “missing” indicators to the network. This way, the model can learn possible patterns from missing data.
if transform.missing_cols.size == 0:
missing = None
else:
missing = np.isnan(x[:, transform.missing_cols]).astype(np.float32)
When In Doubt, Embed
This next fix is small, but it could benefit our model quite a lot. Previously, we’d been feeding all the numerical features directly into one dense layer. Essentially, all the standardized numerical values were stacked together into one long vector and fed to the first layer of the network. The network certainly got all of the data, but it also mixed everything all at once. This forced the first layer to understand what each feature meant and how features should interact. That can be a lot to ask on a first encounter!
The new iteration changes this by using feature-wise numerical embeddings. Before being combined with other features, each numerical feature is first given its own small learned representation. So instead of treating each feature as a single number, the model passes that one value through a tiny transformation to convert it into a short vector. This happens independently for each feature. For example, flux gets its own embedding, spectral curvature gets its own, each energy bin gets its own, and so on. Only after these individual embeddings are they concatenated and passed into the shared model trunk.
While this seems like an extra step, it lets the model learn each feature on its own before combining them. Some features have nonlinear importance, so feature-wise embeddings allow the model to capture that early. This way, the model can learn more meaningful information from preprocessed feature representations.
Aggregation Is Key
The final change I’ve introduced to the model relates to how we obtain the final prediction and uncertainty. This important fix could materially impact the evaluation metrics because they depend on the predictions and uncertainties.
At this point, I must admit that I made a mistake…but don’t we all? While looking back through last week’s code, I noticed that my model had been overwriting predictions on every fold. The setup used fifty global folds for repeated cross-validation so that each source could be held out and predicted multiple times across different repeats. However, the original model would take the latest prediction and overwrite the old one. This is incorrect because repeated cross-validation is specifically designed so that each source gets evaluated multiple times under slightly different training and validation splits. We’re losing a lot of information by throwing away older predictions.
The correct way to do this would be to average all of its repeated held-out predictions, not just the final one. Not only does this provide a more stable prediction, but it also handles uncertainty more accurately. The original model only captured the model uncertainty, but we also need the uncertainty across different repeats. Since each repeat trains on a different outer training set, predictions will also vary across each repeat. Our model needs to be able to capture this uncertainty as well, which is why we need this new fix. The new model uses the Law of Total Variance to explicitly combine both uncertainty sources.
Results
Unfortunately, I don’t have final results to present and make pretty plots with this week because I ran into an error. In the previous section, I talked about aggregating all the held-out predictions. However, since I didn’t use this method for the baseline models, there’s no way to compare the aggregated predictions. Visual Studio Code recognized that before me and angrily spit out this error message at me:
------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[9], line 52
39 rf_aligned, rf_proba, rf_std, rf_counts =
40 align_repeated_summary_to_train_df(
41 train_df=train_df,
42 agg_df=rf_oof,
43 prob_cols=PROB_COLS,
44 model_prefix="rf",
45 )
46 cb_aligned, cb_proba, cb_std, cb_counts =
47 align_repeated_summary_to_train_df(
48 train_df=train_df,
49 agg_df=cb_oof,
50 prob_cols=PROB_COLS,
51 model_prefix="cb",
52 )
---> 53 validate_repeated_counts("RandomForest", rf_counts,
54 expected_counts)
55 validate_repeated_counts("CatBoost", cb_counts,
56 expected_counts)
57 print("=== Data status ===")
Cell In[8], line 148
146 mismatch = np.where(observed != expected)[0]
147 first = mismatch[:10]
---> 148 raise ValueError(
149 f"{model_name}: repeated-CV validation counts disagree
150 for {len(mismatch)} sources. "
151 f"First mismatches at rows {first.tolist()}:
152 observed={observed[first].tolist()}
expected={expected[first].tolist()}"
153 )
ValueError: RandomForest: repeated-CV validation counts disagree
for 2877 sources.
First mismatches at rows [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]:
observed=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
expected=[10, 10, 10, 10, 10, 10, 10, 10, 10, 10]
I haven’t had the time to address this hiccup yet, but I assume that there’s a relatively straightforward fix. I’d probably have to apply the same aggregation pipeline to the baseline models and re-evaluate them. It definitely sounds a little tedious, but I want to aim for as much accuracy and integrity as possible here. I’ll continue working on it and report my results in next week’s post. Until then, stay curious, my fellow reader!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
BNN – Bayesian neural network
SED – spectral energy distribution
References
[1] Gorishiniy et al. “Revisiting Deep Learning Models for Tabular Data.” Conference on Neural Information Processing Systems. 2021. arXiv, doi.org/10.48550/arXiv.2106.11959
Week 11: Ground Control to Major Progress
Hi there, and thanks for stopping by again for the second-to-last week of my senior project! This week, I can confidently say that we’ve made some major progress. If you recall, we unfortunately ended with an error in my code last week. But as they all say, you need a little bit of darkness to see the light. That certainly is true, for this week, the light at the end of the tunnel has become much brighter. I won’t hold the good news from you any longer, so let’s get right into it!
One Error Later
First things first, we have to fix the error from last week. The rather long error code looked like this:
------------------------------------------------------------
ValueError Traceback (most recent call last)
Cell In[9], line 52
39 rf_aligned, rf_proba, rf_std, rf_counts =
40 align_repeated_summary_to_train_df(
41 train_df=train_df,
42 agg_df=rf_oof,
43 prob_cols=PROB_COLS,
44 model_prefix="rf",
45 )
46 cb_aligned, cb_proba, cb_std, cb_counts =
47 align_repeated_summary_to_train_df(
48 train_df=train_df,
49 agg_df=cb_oof,
50 prob_cols=PROB_COLS,
51 model_prefix="cb",
52 )
---> 53 validate_repeated_counts("RandomForest", rf_counts,
54 expected_counts)
55 validate_repeated_counts("CatBoost", cb_counts,
56 expected_counts)
57 print("=== Data status ===")
Cell In[8], line 148
146 mismatch = np.where(observed != expected)[0]
147 first = mismatch[:10]
---> 148 raise ValueError(
149 f"{model_name}: repeated-CV validation counts disagree
150 for {len(mismatch)} sources. "
151 f"First mismatches at rows {first.tolist()}:
152 observed={observed[first].tolist()}
expected={expected[first].tolist()}"
153 )
ValueError: RandomForest: repeated-CV validation counts disagree
for 2877 sources.
First mismatches at rows [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]:
observed=[1, 1, 1, 1, 1, 1, 1, 1, 1, 1]
expected=[10, 10, 10, 10, 10, 10, 10, 10, 10, 10]
The final few lines indicate that the error arises from a mismatch in the repeated cross-validation (CV) counts. This is because last week, we fixed how we calculated the held-out predictions: rather than overwriting predictions in every repeat, we aggregated all predictions and averaged at the end. However, we didn’t do that originally for our baseline models, so they don’t follow the corrected scheme. This means they have a single prediction value (observed=1 in the error message) rather than one prediction for every repeat (expected=10 in the error message).
At this point, I really only had two solutions. I could either completely retrain both baseline models with the corrected prediction calculation, or I could simply change how I compared the models. I didn’t want to take any shortcuts, so I decided to retrain the baseline models, as tedious as it seemed. After retraining and saving the new aggregated predictions, I successfully loaded the baseline results into my code and proceeded to train the improved Bayesian neural network (BNN).
A Not-So-Starbucks Refresher
Since I’ve already described all of the fixes we’re making in this new BNN iteration in last week’s post, I won’t go into detail this time. I’ll just provide a brief list of all the improvements for a quick recap.
Missing data: Previously, we median-imputed Block B but didn’t expose missing data or values. For tabular data, the fact that a value is missing can carry information. Here we keep the imputed/scaled values and pass binary “missingness” indicators to the network.
Feature-wise embeddings: Instead of feeding standardized scalars directly into one dense layer, we’re embedding each numerical feature separately and only then mixing them in the shared trunk.
Correct repeated-CV aggregation: The earlier iteration overwrote predictions fold by fold. Here, we fix it by aggregating and averaging both the repeated-CV predictive mean and the uncertainty.
Model reconstruction: Of course, the biggest change we made last week was reconstructing the entire model. Rather than the computationally expensive, three-model architecture, we’ve produced a shared-trunk hierarchical BNN. This only requires one model but still includes a robust hierarchical strategy with Bayesian branches in the final layers.
Drumroll for the Results
After making the aforementioned changes, I’ve arrived at some very promising results! As always, let’s compare the evaluation metrics of the new BNN to those of the baseline models and the old BNN. The table below shows the six metrics we’ve become familiar with over the past several weeks:
Accuracy | F1 | ROC-AUC | NLL | Brier | ECE
Retrained RF | 0.8495 | 0.8419 | 0.9522 | 0.4276 | 0.0442 | 0.0475
Retrained CB | 0.8442 | 0.8429 | 0.9507 | 0.4273 | 0.0446 | 0.0126
Old BNN | 0.7810 | 0.7762 | 0.9328 | 0.5450 | 0.0570 | 0.0285
New BNN | 0.8443 | 0.8466 | 0.9541 | 0.4485 | 0.0459 | 0.0184
From this table, we can clearly see that the repaired BNN has made some considerable progress! It performs competitively against both baseline models across all metrics, even outperforming them on a couple. Compared to the old BNN, the results are pretty dramatic. Both the accuracy and F1 score are now around 0.85. This shows that the new model isn’t only making more predictions overall, but also balancing precision and recall much more effectively across classes.
On the uncertainty side, the new BNN has also substantially reduced its negative log-likelihood (NLL) and Brier score relative to the old version, indicating that its probability estimates are far more trustworthy and better calibrated. Its expected calibration error (ECE) has dropped to around 0.02, putting it much closer to the highly calibrated CatBoost model. Overall, the new BNN is finally behaving like a serious competitor!
Pretty Pictures
There’s no better way to put results into perspective than with some plots. This week, I’ve prepared both calibration and uncertainty plots for you. The calibration plots, shown below, compare the performances of the retrained baseline models with the new BNN. They paint a pretty encouraging picture (ha, see what I did there).
Across most classes, the new shared-trunk BNN tracks the ideal diagonal line very closely, meaning its predicted probabilities generally match reality quite well. In the blazar-like classes (bll_like and fsrq_like), all three models perform strongly, but the BNN and CatBoost tend to stay slightly closer to the diagonal at high probabilities. The young_psr_like class also looks noticeably improved for the BNN, with an ECE competitive with the other models. The one trouble spot remains other_galactic; all models continue to struggle there due to the tiny sample size and severe class diversity. Still, compared to the old BNN from previous weeks, these curves are a major step forward.
The next few plots we have capture the models’ uncertainty measures. Looking at the very left panel, we can see that the new BNN captures uncertainty much better than before. Compared to previous weeks, the uncertainty is more strictly and negatively correlated with the max probability. In other words, the new model is less uncertain for higher predicted probabilities, which is exactly what we want. Using the color-coding legend, we can also see that there are fewer incorrect predictions at higher probabilities than before.
Now, let's turn our attention towards the two plots on the right. Both plots compare the uncertainties of the three models. From the center plot, we can see that both baseline models’ uncertainties spike near zero, while the new BNN’s uncertainty is more spread out. The plot on the right confirms this, revealing that the BNN’s average uncertainty is higher across all classes. While this seems discouraging, it’s not completely unexpected or something to worry about too much. Since the baseline models don’t have uncertainty integrated into their model architectures, they only capture the variability between cross-validation repeats. On the other hand, the BNN incorporates the posterior predictive uncertainty as well. The BNN is therefore measuring something much broader, so a higher uncertainty doesn’t automatically mean it’s worse. That being said, it is definitely something we should look into next week!
Looking Ahead
Now that we’ve got more optimistic results, we can focus on refining and confirming the small details. In the table a couple of sections above, notice how I put the old-BNN numbers in red. I did that on purpose because those numbers aren’t entirely honest. That’s not to say they’re “wrong” in any way, but they just don’t tell us the full story. Unlike the other three models, the old BNN doesn’t use the correct prediction aggregation; it’s still using the old, overwritten-prediction scheme. For reliable comparisons, we should use the same aggregation strategy on the old BNN while keeping everything else the same.
Next week, I’ll focus on two things. The first will be to apply each change we discussed last week, one at a time, to the old BNN. This way, we’ll be able to see which changes create the most substantial improvements. Even though the new BNN performs quite well, it’s worth checking each fix to see if any of them are nonessential to the model’s enhancement. That would 1) give us a better idea of which architectural and training choices are actually driving the performance gains, and 2) possibly save us some computational complexity.
The second key area, which I mentioned when discussing the plots, will be to investigate the new BNN’s uncertainty. I have a couple of ideas in mind, which I’ll save for next week’s blog post. Stay tuned for the final week, where we’ll see how our ultimate model iteration performs. Until then, stay curious!
Acronyms
Here’s a list of all the acronyms used throughout this paper. I ordered them in terms of first appearance:
BNN – Bayesian neural network
CV – cross-validation
NLL – negative log-likelihood
ECE – expected calibration error
Week 12: Messages from the Final Frontier
Hello, my fellow reader! It is with great joy that I welcome you to the final blog post of my senior project. I spent most of this past week tidying up my results for my paper and presentation. I’ve made a few small but important changes since I last spoke with you, so let me break them down and show you the closing results of this project!
John: One Fix At A Time
At the end of last week’s post, I mentioned that I would be seeing what happens if I applied changes to my model one at a time. Rather than running all of the changes to my original BNN at once, I’ve created two new iterations. The first only fixes the prediction aggregation, and the second includes all the minor improvements except the major architectural change. By splitting up the changes a little, I’ve been able to better understand which ones create the most improvement for my model.
Prediction Aggregation
With just the aggregation fix, my model’s evaluation metrics jump up significantly. This is as expected because, originally, the predictions were overwriting each other during training with every consecutive repeat. This made the repeated cross-validation process effectively equivalent to running only a single fold, which is not very robust for a model! Once all ten repeats are properly accumulated and averaged together, the model finally gets a much more stable and reliable estimate of its performance.
Architecture Comparison
Next, I’ve applied all of the changes I made last week, besides the architectural change. In other words, I’m using the three-model, hierarchical BNN, but incorporating all of the training and evaluation fixes we discussed last week. This essentially tells us whether the new, shared-trunk BNN architecture has any real benefit.
Now that we’ve applied the separate changes, we can compare their effects on our model in the table below. The first column shows the evaluation metrics for the original hierarchical BNN without any changes. The next column shows how performance changes when the predictions are aggregated, while keeping the same architecture. The third column then reflects the impact of all the other non-architectural modifications applied to the hierarchical BNN. Finally, the fourth column compares the shared-trunk hierarchical BNN with all of the changes that we made last week.
OG BNN | OG BNN Aggregated Predictions | OG BNN All Other Changes | Shared-Trunk BNN
Accuracy | 0.7810 | 0.8120 | 0.8186 | 0.8384
F1 | 0.7762 | 0.8068 | 0.8071 | 0.8413
ROC-AUC | 0.9328 | 0.9406 | 0.9502 | 0.9530
NLL | 0.5450 | 0.5289 | 0.4882 | 0.4584
Brier | 0.0570 | 0.0559 | 0.0515 | 0.0470
ECE | 0.0285 | 0.0523 | 0.0412 | 0.0343
Looking at the table, we can see that the original BNN performs worst on nearly all of the metrics. Each new iteration makes small improvements over the last, which is exactly what we want to see. In the end, the shared-trunk BNN outperforms all three hierarchical variants, which tells us that: 1) the seemingly minor training and evaluation fixes we made actually did make a difference, and 2) the general learned representation of the shared-trunk BNN is superior to the isolated pathway approach of the hierarchical BNN.
Paul: Uncertain About Uncertainty?
Moving on to our next update, I’ve decomposed the shared-trunk BNN’s uncertainty. I’m doing this so that I can better compare the uncertainty to that of the baseline models. Because the tree-based baseline models are non-Bayesian, they don’t have that intrinsic uncertainty measure that comes from the posterior in Monte Carlo sampling. That means the BNN has an extra component of uncertainty, making the comparison unfair. So instead of using a single scalar to quantify the BNN’s uncertainty, I’ve split it up into two components.
Between-Repeat
The first is the between-repeat variance, which comes from how the data is randomly split for every cross-validation repeat. In every repeat, the model learns from a slightly different subset of training data, so there’s going to be some variation in its predictions. If we take the standard deviation of the predictions across all repeats, we can then derive an uncertainty quantification from this repeated cross-validation strategy. Since we use the same cross-validation scheme for the tree-based baselines, this uncertainty component is the one that’s actually comparable across all models.
Posterior Predictive
The second component is the built-in BNN uncertainty, which comes from the stochastic nature of the Bayesian heads in the shared-trunk hierarchical architecture. During inference time, the model is run 200 times with different weight samples drawn from its learned posterior, producing a distribution of predictions for each source. The variance across these MC draws captures how uncertain the BNN model’s own weights are, independent of how the data was split. This is unique to the BNN and has no equivalent in the tree-based models, which makes it harder to compare.
Now that we’ve split the uncertainty, we can more reliably compare it across the three models. In the first plot below, I’ve shown the average uncertainties for the different models. The left panel is the between-repeat variance that’s directly comparable across each model. We can immediately see that, even with the decomposed uncertainty, the BNN is still more honest. Both the Random Forest and CatBoost models underestimate their uncertainties. Interestingly, the posterior uncertainty tends to be slightly lower than the between-repeat uncertainty, but the two track each other pretty closely in terms of which classes the model is most uncertain about. All of the models seem relatively confident about the dominant bll_like and fsrq_like classes, but they remain more uncertain about smaller classes like msp_like and other_galactic.
In this next plot, we can see the per-source uncertainties of each tree-based model plotted against those of the shared-trunk BNN. Essentially, it shows us the correlation between the baselines and the BNN. The uncertainty points for both the Random Forest and CatBoost models lie below the diagonal. This indicates that they’re overconfident in their predictions. Despite this, the correlations are reasonably strong, suggesting that all three models largely agree on which sources are hard to classify. The coloring of the points further reinforces this, as the incorrectly classified sources (shown in red) tend to sit higher and further to the right. This means that the higher BNN uncertainty is a reliable signal that a source is ambiguous. Overall, the BNN seems to report more uncertainty in the right places than the tree-based models.
George: Looking Ahead
If you’ve been along over these past twelve weeks, you’ll know that I always like to end every post with a short section that previews the following week’s tasks. Even though this is the last week of the senior projects, I still wanted to keep this practice. Except this time, it won’t be merely a look ahead into the coming week. Rather, it'll be a discussion of the future work that could extend this research project.
Given the limited timeframe of this program, I didn’t get the chance to fully accomplish everything I had hoped for. Whether it’s taking completely new directions or making ongoing refinements to my model, there are a lot of ideas I have in mind that would take this project even further.
To improve my model, I’d like to experiment with various data augmentations. I’d first like to incorporate more data, whether that’s different types of data or data from different observatories. My project only used spectral and time variability data, but in the future, I’d definitely like to see how my model performs with additional types or instances of data. Along with that, I’d like to perform some data augmentation on my existing data. We could then see if that could account for the severe class imbalance. Previous papers have shown that this might not be the best solution[1], but you’ll never know until you try for yourself, right?
Once I’ve stabilized my BNN architecture, I’ll consider training more probabilistic models. I’m currently thinking about Gaussian Processes and Mixture Density Networks, but I’ll definitely look into more algorithms. This would be a good way to see how different probabilistic methods quantify uncertainty, and whether more complex architectures are better for this classification problem. The results of this study would be significant in helping researchers determine which probabilistic approach is best suited for classifying ambiguous gamma-ray sources.
Finally, but perhaps most importantly, I’d like to actually apply my model(s) to classify unassociated sources in the Fermi-LAT catalog. If you recall, the vast majority of sources in the Fermi-LAT 4FGL-DR4 catalog[2] remain unassociated. That was exactly what my project sought to tackle in the first place. Unfortunately, I didn’t get the time to accomplish that goal within these limited twelve weeks, but that simply leaves more exciting ground to cover in the future! Doing so would meaningfully advance our understanding of the high-energy gamma-ray sky, and potentially reveal new source populations in the process.
Ringo: Reflections
Before I bid you farewell for good (I know, I know, you’re getting tired of me), I wanted to take some time to reflect on this experience and to give my proper thanks to the people and resources that have supported me along the way. If you haven’t already noticed, the headings of each section in this post are the names of each Beatles member. This choice was inspired by that one scene in Project Hail Mary, where Dr. Grace sends his findings back to Earth using four space probes, each named after a different Beatle. Just like those probes carrying their discoveries across the void of space, each section of this post carries a piece of the work and thinking that went into this project over the past twelve weeks. It has been an incredible learning journey, with so many successes and challenges to reflect on.
These past several weeks have simultaneously felt like the longest and shortest period of my life. It felt short because I had so much fun devising, planning, and executing my project, but it also felt long because of all the requirements, deadlines, and lengthy blog posts. I’ve literally never written so much and so frequently in my life, but I am not complaining about that at all. In fact, I am so grateful to have had the opportunity to lean into my love for writing while documenting my project for you all to see.
Lesson Learned
Aside from the writing, I thoroughly enjoyed the continuous process of learning, succeeding, failing, and learning again. This project was in no way a linear climb, and I faced numerous challenges throughout. However, I think that might be the most valuable outcome of this project. Looking back on these twelve weeks, I’ve noticed both academic and personal growth. Whether it was searching (and asking Claude) to understand tedious equations or scrutinizing my code to find that one small error, every obstacle I pushed through made me a little bit better at this than I was before. I’ve definitely become a more knowledgeable and resilient researcher in the process.
For the juniors who are hopefully reading this, I advise you to pursue a passion project throughout your senior year. Even if you’re not asking for all of the structuredness and supervision of the senior project, you'll gain so much from dedicating yourself to something you genuinely care about. In my very personal and non-biased opinion, I think the senior project is a great option. While it can be demanding at times, it’s an incredibly rewarding program that gives you dedicated time and resources to see a project through from start to finish. It’s a rare opportunity for you, a literal high schooler, to go out, meet professionals in your field, expand your network, and contribute to something you truly love. Plus, you’ll get to graduate with honors!
Acknowledgements
Yes, I’m very excited to graduate with honors, but I’m more proud of all the wonderful relationships I built with the people involved in this process. This project was so much less about the results and so much more about the people who made it memorable. It wouldn’t have been nearly possible without the support of so many brilliant and generous individuals. I’d like to take this time to thank each and every one of them.
First off, I’d like to thank my parents for always having my back and making the 2-hour round trip to SFSU every week. Without their benevolent Ubering services, I wouldn’t have been able to meet my external advisor, Dr. Oscar Macías, who has been an incredible mentor throughout this entire process. I’m so grateful that I was able to work with a professor who was so enthusiastic and knowledgeable, and always willing to point me in the right direction.
I’d also like to extend my thanks to my internal advisors, Mrs. Bhattacharya and Mr. Schindler, who were just as instrumental in getting me here. Between keeping me on track with deadlines, talking through technical difficulties, and offering emotional encouragement, their support was something I could not have done without. Of course, I’d like to acknowledge the school for providing the opportunity that made this project possible. And finally, I’d like to thank you, my dedicated reader, for bearing with me for this long.
With that, I’ll bid my final farewell. As TARS from Interstellar perfectly puts it, I’ll “see you on the other side.” It’s been an absolute pleasure, and until next time, don’t forget to look up!
Where To Find Me
If you were unable to make the final presentations, you can find the presentation recording linked below. Along with that, I've also provided the links to my presentation slides and GitHub repository. Feel free to check them out!
Presentation recording: https://drive.google.com/file/d/1e4NTtqUqvNhJtR-2jPlXgJjPDoU85rLP/view?usp=sharing
Presentation slides: https://canva.link/8x1hsvb62ehgut5
GitHub repository: https://github.com/chloelui/Fermi-Pulsar-BNN
References
[1] Butter et al. “Classification of Fermi-LAT blazars with Bayesian neural networks.” Journal of Cosmology and Astroparticle Physics, 2022. IOPscience, 10.1088/1475-7516/2022/04/023.
[2] NASA Fermi-LAT Collaboration. (2024). LAT 14-year Source Catalog (4FGL-DR4). Version 4, Fermi Science Support Center, 24 Jul. 2024, fermi.gsfc.nasa.gov/ssc/data/access/lat/14yr_catalog/. [Dataset].
[3] Boer, Poppe de. “The Beatles in Treslong.” Wikimedia Commons, 5 Jun. 1964, commons.wikimedia.org/wiki/File:The_Beatles_in_Treslong.jpg. Public Domain.