Rock, carbon or metal: what are asteroids made of?¶
EPS 88 · PyEarth · final project. Open your own copy on DataHub: click here.
Nobody has visited (16) Psyche, an asteroid that may be made largely of metal, so what it is made of was read from afar: how much sunlight it reflects and what colour it is. In this project you teach a computer to read the same two things and sort asteroids into rocky, carbon-rich, and the group that holds the metal worlds. Then it sorts 368 classified asteroids whose classes you never see. The leaderboard sees them, and it ranks everyone's entries.
This notebook is built like a weekly notebook: a few questions lead it, each section answers one, and every place you work opens with a pencil icon ✏️ and the words Your turn. The difference is how much you write. Each Your turn says why the step matters and gives a numbered plan in words; the cell under it is empty, and you write all of the code, looking back at the weekly notebook it names. A self-check cell follows. Run it: it tells you whether your numbers are right, as in every week.
The last Your turn is your leaderboard entry: it saves the file submission.csv. Run every cell so your answers and figures are saved in the file, then download this notebook itself — the .ipynb file — and upload it to Gradescope together with submission.csv. In JupyterLab: File ▸ Download, or right-click the file in the left-hand panel and choose Download. The leaderboard scores your file on answers you never see and ranks everyone's entries.
Two habits from the first minute. A cell runs when you press Shift+Enter, and the notebook remembers everything it has already run — so when something breaks and you cannot see why, Kernel → Restart Kernel and Run All Cells throws the memory away and rebuilds it from the top. That is never the wrong thing to do.
What your project is about: rock, carbon or metal?¶
Seven steps before any code: the problem, why it matters, the data, and what is known.
1. A spacecraft is flying to an asteroid that may be the metal core of a planet's building block¶
NASA's Psyche spacecraft launched on 13 October 2023 and is due to begin exploring the asteroid (16) Psyche by August 2029. Psyche orbits the Sun between Mars and Jupiter and is about 280 km across at its widest. Scientists think it holds a lot of metal and may be part of the core of a planetesimal, one of the small bodies the planets were built from (NASA). The mission is led by Lindy Elkins-Tanton, a professor in Berkeley's Department of Earth and Planetary Science.
How Psyche might look if metal (grey-blue) and rock (brown) sit in large patches. No spacecraft has photographed it yet. Illustration: NASA/JPL-Caltech/ASU, PIA24471.
2. Why it matters: nobody can visit a planet's core¶
Earth's core lies under the whole rocky mantle, out of reach; Psyche may show us a core. NASA says Psyche is most likely the survivor of several violent hit-and-run collisions from when the Solar System formed, so it may tell us how the cores of Earth and the other rocky planets came to be. Asteroids matter for a second reason: they are leftover building blocks of the planets, and the meteorites that land on Earth come from them. Their colours have already tied several kinds of meteorite to the asteroids they came from (Reddy and others, 2015, arXiv:1502.05008).
Earth's layers, drawn to scale on the right: the core is the grey centre. Diagram: Surachit, Wikimedia Commons (CC BY-SA 3.0).
3. An asteroid is read from its colour¶
An asteroid is a point of light, so astronomers read its make-up from the colours of that light. Each mineral absorbs some colours of sunlight more than others. The rock minerals olivine and pyroxene absorb near 1,000 nm, in the near infrared just beyond what our eyes see (Reddy and others, 2015, arXiv:1502.05008). So astronomers measure an asteroid's reflectance spectrum: the light it reflects at each wavelength, divided by the light it reflects in green (550 nm). They sort spectra into classes by shape, and group the classes into three complexes. This notebook calls each complex a kind:
- S (stony): rises steeply toward the red and dips near 1,000 nm, the dip of olivine and pyroxene. Gaia's colours (step 5) stop at about 1,000 nm, so they see only the start of this dip, where an S spectrum begins to fall after about 800 nm;
- C (carbon-rich): nearly flat;
- X: rises gently with no dip. Psyche is an X, but so are asteroids that only share its colour.
Spectrum shapes of 12 of the 24 classes of the Bus–DeMeo system (2009). Across: wavelength from 0.45 to 2.45 µm (the "1" marks 1,000 nm). Up: reflectance, set to 1 at 0.55 µm. Compare S, C and X. Diagram: Wiki Remi, Wikimedia Commons (CC BY-SA 4.0).
4. Why a model: telescopes have classified only a few asteroids¶
A spectrum good enough to classify an asteroid needs a large telescope and a bright asteroid, so only a few have one. Your training file holds 138,094 asteroids with a measured albedo. Gaia measured the colours of 40,642 of them, but only 1,534 have a class from telescope observations, and 861 have both a kind and a full Gaia colour spectrum. A model can learn from those few how albedo, orbit and colour go with the kind, and then name a kind for the tens of thousands that no telescope has classified.
Counted from the training file of this project. The axis is logarithmic: each step is ten times more asteroids.
5. The data: two space telescopes, joined on the asteroid's number¶
One row of the data is one asteroid, and 97% of the rows are in the main belt, between Mars and Jupiter. NASA/JPL's Small-Body Database gives each asteroid's orbit, size, albedo and class letters. Psyche, for example, orbits 2.93 AU from the Sun (1 AU is the Earth's distance), is 222 km across on average (it is not round, so its widest is longer, step 1), and has an albedo of 0.12.
The inner Solar System: the main belt is the white ring between the orbits of Mars and Jupiter. Image: Mdf, Wikimedia Commons (public domain).
- Albedo, the share of sunlight a surface reflects (0.1 means a tenth), comes from infrared telescopes. It is measured head-on, with the Sun behind the observer, and compared with a plain white disc; a very bright surface can reflect as much as that disc or more, so a few rows of the file reach 1. NASA's WISE measured the heat glow of more than 100,000 main-belt asteroids, which gives each one's size; size and visible brightness together give the albedo (Masiero and others, 2011).
- Colour comes from ESA's Gaia spacecraft, built to map the stars of the Milky Way, which also watches asteroids. Its third data release gives reflectance spectra of 60,518 asteroids in 16 bands from 374 to 1,034 nm, observed from 2014 to 2017: each is the asteroid's light divided by the light of stars like the Sun (Galluccio and others, 2023).
- Classes are letters from two systems built from telescope observations: Tholen's (1984), from colours and albedo, and SMASSII (Bus and Binzel, 2002), from spectra.
Left: NASA's WISE in orbit around Earth (NASA/JPL-Caltech, PIA17254, public domain). Right: ESA's Gaia (ESA, CC BY-SA 3.0 IGO). Both artist's impressions, via Wikimedia Commons.
6. What earlier research found¶
For more than forty years astronomers have found that the kinds are not spread evenly through the belt.
- Gradie and Tedesco (1982) sorted asteroids into colour classes and found that the share of each class changes steadily with distance from the Sun, from 1.8 to 5.2 AU. They read this as asteroids having formed near where they orbit now. doi:10.1126/science.216.4553.1405
- Masiero and others (2011) measured the size and albedo of more than 100,000 main-belt asteroids with WISE. Albedo has two humps, dark and bright, and the average albedo changes from the inner to the outer belt. doi:10.1088/0004-637X/741/2/68
- DeMeo and Carry (2014) mapped the belt with many smaller asteroids and found it more mixed, by size and distance, than had been known: a sign that the belt was stirred, for example as the planets moved. doi:10.1038/nature12908
- Delbo and others (2026) trained a classifier on Gaia spectra of asteroids with known classes and named a class for every asteroid in the Gaia release. S was sorted well; X mixed with its neighbours, as Gaia's short range of wavelengths leads one to expect. doi:10.3389/fspas.2026.1774478
The share of each class against distance from the Sun, redrawn from Gradie and Tedesco (1982). S (red) peaks in the inner belt and C (blue) farther out; the grey band is the main belt. Chart: Praemonitus, Wikimedia Commons (CC BY-SA 4.0).
7. What is still open, and what you will try¶
Whether an X asteroid is metal is still open. The M class (part of X) was first named for asteroids thought to be the metal cores of broken planetesimals. Psyche, the best-known M-type, was once thought to be almost all metal. New measurements of its density, its radar echo and its spectrum point instead to a mix of metal and rock: Elkins-Tanton and others (2020) estimated a density of 3,400 to 4,100 kg/m³ and 30 to 60% metal by volume, and wrote that final answers may wait until the spacecraft arrives (doi:10.1029/2019JE006296). A newer classification found that albedo separates X asteroids whose spectra look alike (Mahlke and others, 2022, doi:10.1051/0004-6361/202243587).
What you will try: build a classifier that sorts S, C and X from albedo, orbit and Gaia colour; test it on asteroids it has never seen; then ask whether it can be trusted on the small asteroids that no telescope has classified.
NASA's Psyche spacecraft and its two cross-shaped solar arrays. Illustration: NASA (public domain), via Wikimedia Commons.
What you'll be able to do¶
The science. Say why rocky asteroids look bright and red, carbon-rich ones dark and grey, and the metal-bearing group in between; why the rocky ones crowd the inner asteroid belt; and why colour, not orbit, is what tells the kinds apart.
The code. You write every line: a histogram for each group, a function for the macro F1 score,
a rule with two cuts, RandomForestClassifier, a stratified split, a bootstrap, an SVM with
scaled columns, two columns made from a colour spectrum, and a submission.csv file.
The four questions of the project:
- How bright are rocky, carbon-rich and metal-bearing asteroids?
- Can brightness and orbit sort the three kinds?
- How well does the classifier do on asteroids it has never seen?
- Does a different model, or more information, sort them better?
What you hand in, and when¶
One Gradescope assignment, "Final project: asteroids", takes everything.
- Individual work. Start now; Gradescope takes uploads from Wed 25 Nov, and the project is due Wed 16 Dec. Mon 7 Dec (RRR week) has no class, but you are welcome to come to the room in class time and work on your notebook.
- Every upload carries two files:
submission.csv(your leaderboard entry, made by the last Your turn) and this notebook. Gradescope cannot check this for you, so follow these steps:- In DataHub, save the notebook (File ▸ Save).
submission.csvis in the same folder as the notebook, in the left-hand panel. - Download both: right-click each file in the left-hand panel and choose Download.
- In Gradescope, open Final project: asteroids, click Upload Submission, and select both files (hold Ctrl, or Cmd on a Mac, to select the second), or drag both into the window. Type a name for the leaderboard: your own name or a nickname. Your classmates will see it next to your score.
- Check that Gradescope lists both files, then upload. If the result says a file is missing, upload both again.
- In DataHub, save the notebook (File ▸ Save).
- Up to 10 uploads that change
submission.csvare scored on the leaderboard. An upload that changes only the notebook is free. - Everything is due Wed 16 Dec at 23:59 Pacific time: the notebook and your last entry.
- Before every upload: Kernel → Restart Kernel and Run All Cells must run from top to bottom without an error.
How it is graded (100 points). The Your turns are graded like a week's. The leaderboard gives two grades: one for a valid entry that beats always guessing the most common kind, and one for your rank, on the final score from the hidden part of the board. Choose your entry by your own fair test, not by the public score.
| What | Points | Graded by | Full marks when |
|---|---|---|---|
| Your turns 1 to 10 | 80 | Gradescope, when you upload | each Your turn by its self-check (Your turn 1, a figure, for being done and run), and the notebook runs from top to bottom without an error |
| Leaderboard entry | 10 | Gradescope, pass or fail | your last submission.csv is valid and beats always guessing the most common kind |
| Leaderboard rank | 10 | after the board closes | your final score, on the hidden part, ranked among the students on this project: top quarter 10, second quarter 8, third quarter 6, last quarter 4; a tie takes the higher band; no valid entry 0 |
Setup¶
Run the next cells once. All their code is given. One row is one asteroid (step 5 at the top). Three files ship with the course:
trackT14_asteroid_catalogue.csv.gz: every numbered asteroid with a measured albedo, your training file, frozen so everyone trains on the same file. The 368 leaderboard asteroids have been taken out of it;trackT14_leaderboard_inputs.csv: the leaderboard file, those 368 asteroids with their orbit, size, albedo and colour but not their class, and no name;trackT14_submission_template.csv: the empty submission template, one row per leaderboard asteroid.
Use only these files. The online catalogues hold the hidden classes, and the leaderboard flags a notebook that reads them.
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.metrics import f1_score
from sklearn.model_selection import train_test_split
from sklearn.pipeline import make_pipeline
from sklearn.preprocessing import StandardScaler
from sklearn.svm import SVC
import pandas as pd
import matplotlib.pyplot as plt
# house style, set once, so every plot cell below holds only what matters
plt.rcParams.update({"figure.figsize": (7, 4.5), "figure.dpi": 110,
"axes.grid": True, "grid.alpha": 0.3, "axes.axisbelow": True})
CACHE = "https://raw.githubusercontent.com/AI4EPS/EPS88_PyEarth/main/data"
# larger text for every figure, and one colour per kind, used in every figure below
plt.rcParams.update({"font.size": 11, "axes.titlesize": 12, "legend.fontsize": 10})
KIND_COLOURS = {"S": "tab:red", "C": "navy", "X": "tab:green"}
asteroids = pd.read_csv(CACHE + "/trackT14_asteroid_catalogue.csv.gz", dtype={"spec_T": str, "spec_B": str})
leaderboard = pd.read_csv(CACHE + "/trackT14_leaderboard_inputs.csv")
template = pd.read_csv(CACHE + "/trackT14_submission_template.csv")
print("training file:", asteroids.shape)
print("leaderboard file:", leaderboard.shape)
What the columns hold.
a,e,i: the orbit.ais the average distance from the Sun in AU (the main belt lies between about 2.1 and 3.3 AU);esays how stretched the orbit is (0 is a circle);iis its tilt in degrees.diameterin km, andalbedo.R374toR1034: Gaia's reflectance in 16 bands, named by wavelength in nm, each divided by the reflectance at 550 nm (green). SoR550is always 1, a value above 1 means "redder than green", and a blank is a band Gaia flagged as unreliable.spec_Tandspec_B: the class letters of the Tholen and SMASSII (Bus) systems. Most asteroids have neither.
The next cell turns the class letters into the three kinds S, C and X of step 3 at the top. It reads the Bus letter first and the Tholen letter if Bus has none, the leaderboard's rule. A letter that names none of the three (D, T, an uncertain class with a colon) gives no kind. Tholen's dark P class counts as X here, as the leaderboard counts it. That is a choice: the newer classification of Mahlke and others (2022, step 7) puts P with C instead. In the data, 13 of the 194 X asteroids are X only through Tholen's P. Then the cell keeps the asteroids with a kind and all 14 bands from 418 to 990 nm, because every leaderboard asteroid has those.
# the letters of each kind, in each of the two systems
BUS = {"S": ["S", "Sq", "Sl", "Sa", "Sk", "Sr", "Q", "A", "K", "L", "Ld", "R", "V", "O"],
"C": ["C", "Ch", "Cb", "Cg", "Cgh", "B"],
"X": ["X", "Xc", "Xk", "Xe"]}
THOLEN = {"S": ["S", "A", "Q", "R", "V"],
"C": ["C", "B", "G", "F"],
"X": ["X", "M", "E", "P"]}
def kind_of(letter, system):
"""S, C or X for one class letter of one system, or None if it names none of them."""
for kind, letters in system.items():
if letter in letters:
return kind
return None
# Bus first; Tholen only where Bus gives no kind
kinds = []
for bus, tholen in zip(asteroids["spec_B"], asteroids["spec_T"]):
kind = kind_of(bus, BUS)
if kind is None:
kind = kind_of(tholen, THOLEN)
kinds.append(kind)
asteroids["complex"] = kinds
# the 14 colour bands every leaderboard asteroid has (550 nm is always 1, but kept for simplicity)
bands = ["R418", "R462", "R506", "R550", "R594", "R638", "R682", "R726", "R770", "R814",
"R858", "R902", "R946", "R990"]
labelled = asteroids[asteroids["complex"].notna() & asteroids[bands].notna().all(axis=1)]
labelled = labelled.reset_index(drop=True)
print("asteroids with a kind:", asteroids["complex"].notna().sum())
print("... and all 14 colour bands:", len(labelled))
print(labelled["complex"].value_counts())
1. How bright are rocky, carbon-rich and metal-bearing asteroids?¶
The first thing to read from an asteroid's light is how much sunlight it reflects: its albedo, the share of sunlight a surface reflects (step 5 at the top). If the three kinds are made of different stuff, they may differ in albedo. Three famous asteroids, one of each kind:
# three famous asteroids, one of each kind
famous = asteroids[asteroids["name"].isin(["Vesta", "Hygiea", "Psyche"])]
print(famous[["name", "complex", "diameter", "albedo", "a"]])
Vesta, the stony one, reflects 0.42 of the sunlight; Hygiea, carbon-rich, 0.07; Psyche, an X, 0.12. Three asteroids are a hint, not evidence. The training file holds 861 asteroids with a kind and a full colour spectrum. Do the three kinds differ in albedo across all of them?
✏️ Your turn 1¶
How is albedo spread within each kind?
Why: before any model, look at the data. If the three kinds sit at different albedos, albedo alone can tell them apart; where their histograms overlap, it cannot.
Look back at: histograms, from Earth's elevation has two peaks (03_two_peaks.ipynb), and the
loop from Which of these worlds could have liquid water? (02_liquid_water.ipynb).
- Make one set of bins for all three histograms: from 0 to 0.6 in steps of 0.02.
- With a loop over the three kinds
"S","C"and"X", take the albedo of the rows oflabelledof that kind and draw its histogram on the shared bins, half transparent (alpha=0.4), in the kind's colour fromKIND_COLOURS, with a label that says the kind, how many asteroids it holds and its median albedo. - Draw the median of each kind as a dashed vertical line in the same colour (
plt.axvline). - Label both axes, add a legend, and give the plot a title with the number of asteroids.
This Your turn has no self-check: a figure you can see needs none. A good figure lets a reader see at a glance where each kind sits: three coloured histograms, three dashed medians, and labels that add up to 861.
# ← your answer here: write your code here
The three kinds sit at three albedos, and they overlap. The middle (median) albedo is 0.23 for S, 0.05 for C and 0.09 for X. S and C barely overlap. X does not have a place of its own: many X asteroids sit with C at the dark end, and the rest spread across the gap and into S.
Think. If you had to name one kind for every asteroid, knowing nothing about it, which would you name?
Scoring a guess with three answers¶
The F1 score you met in the explosion week works on one kind at a time. There it asked two questions of one label, "explosion":
Precision / recall: Of the ones you flagged, how many were right? Of the real ones, how many did you catch?
F1 is the balance of the two. Here there are three kinds. We score each kind the same way, that kind against the other two, and average the three F1 scores. The average is the macro F1. The leaderboard ranks you by it.
Why the average and not the share right? The kinds are not equally common. A rule that never says X can still be right most of the time, and an average of the three F1 scores gives such a rule 0 for X. Macro F1 runs from 0 to 1, and 1 means every asteroid named correctly.
A simple case first: six made-up asteroids, their true kinds and a guess. The next cell scores
one kind, S, with f1_score from the explosion week: a True for every asteroid that is S, against
a True for every asteroid the guess calls S.
# six made-up asteroids: their true kind and a guess for each
toy_truth = np.array(["S", "S", "S", "C", "C", "X"])
toy_guess = np.array(["S", "S", "C", "C", "S", "S"])
# the F1 of one kind, S: that kind against the other two
toy_f1_S = f1_score(toy_truth == "S", toy_guess == "S")
print("F1 for S:", round(toy_f1_S, 3))
The guess called four asteroids S and two of them were; it caught two of the three real S. Precision 2/4 and recall 2/3 give F1 0.571 for S. The same for C gives 0.500, and X, never guessed, gets 0. The macro F1 of the guess is the average of the three.
✏️ Your turn 2¶
How good is the simplest guess: the same kind for every asteroid?
Why: a classifier is worth something only if it beats the dumbest rule, so measure that rule first.
Baseline: Write the dumbest rule you can, first. Any model that cannot beat it is decoration.
Look back at: f1_score, from Earthquake or nuclear explosion? (09_explosion_or_earthquake.ipynb),
and the six made-up asteroids just above.
- Write a function
macro_f1(truth, guess)that takes the true kinds and the guessed kinds and returns the average, over"S","C"and"X", of each kind's F1 score. You will use it for every classifier in this notebook. - Print it for
toy_truthandtoy_guess, and beside it whatf1_scoregives withaverage="macro". The two should agree. most_common: the kind that occurs most often inlabelled.baseline_f1: the macro F1 of namingmost_commonfor every asteroid inlabelled. Print both.
Use these names, because the self-check looks for them: macro_f1, most_common,
baseline_f1.
# ← your answer here: write your code here
# the self-check reads: macro_f1, most_common, baseline_f1
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert abs(macro_f1(np.array(["S", "C", "X", "X"]), np.array(["S", "C", "X", "C"])) - 7 / 9) < 1e-6, \
"macro_f1: the average of the F1 scores of S, C and X, each kind against the other two"
assert most_common == "S", "the kind that occurs most often among the classified asteroids"
assert abs(baseline_f1 - 0.2071) < 0.001, \
"the macro F1 of naming the most common kind for every classified asteroid"
print(f"✓ Your turn 2 — always guessing {most_common} scores a macro F1 of {baseline_f1:.3f}")
Always guessing S scores 0.21. Its F1 for S is 0.62, because S is the largest kind, but C and X get 0, and the average of the three is low. That is the bar: a rule that cannot beat 0.21 has learned nothing about asteroids.
2. Can brightness and orbit sort the three kinds?¶
The histograms suggest a rule with two cuts. Below a low albedo, call it C; above a high one, call it S; in between, X. You choose the two cuts by eye from your histogram of Your turn 1.
✏️ Your turn 3¶
How far does albedo alone get, with two cuts?
Why: this is the simplest rule with any science in it. Every classifier after it has to beat it, or the classifier adds nothing that one number did not already say.
Look back at: np.where, from Earthquake or nuclear explosion? (09_explosion_or_earthquake.ipynb),
and your histogram of Your turn 1.
- Choose two albedo cuts from your histogram:
low_cut, below which most asteroids are C, andhigh_cut, above which most are S. rule_guess: for every asteroid inlabelled,"C"if its albedo is belowlow_cut,"X"if it is fromlow_cutup to belowhigh_cut, and"S"fromhigh_cutup.rule_f1: the rule's macro F1, with yourmacro_f1. Print it besidebaseline_f1.- Try two or three other pairs of cuts and keep the pair with the highest score.
Use these names, because the self-check looks for them: low_cut, high_cut, rule_f1.
# ← your answer here: write your code here
# the self-check reads: low_cut, high_cut, rule_f1
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert 0.02 < low_cut < high_cut < 0.4, "two albedo cuts, the low one below the high one"
again = np.where(labelled["albedo"] < low_cut, "C", np.where(labelled["albedo"] < high_cut, "X", "S"))
assert abs(rule_f1 - f1_score(labelled["complex"], again, average="macro")) < 0.01, \
"rule_f1: the macro F1 of C below the low cut, X between the cuts and S above the high cut"
assert rule_f1 > baseline_f1 + 0.2, \
"the rule should beat always guessing the most common kind by far: are C and S the right way round?"
print(f"✓ Your turn 3 — cuts at {low_cut} and {high_cut}: macro F1 {rule_f1:.3f}")
Two cuts score 0.62, three times the baseline's 0.21. But X is the weak kind: the rule's F1 for X is only 0.27, because the X albedos spread across both cuts. Albedo alone cannot place most X asteroids.
Does where an asteroid orbits say anything more? Earlier research found that the kinds are not mixed evenly (step 6 at the top). The belt itself has gaps, where Jupiter's pull clears asteroids out: the Kirkwood gaps. The next cell splits the belt into three zones at two of them, 2.5 and 2.82 AU, and counts the share of each kind in each zone.
Main-belt asteroids counted by their distance from the Sun a. The five dashed lines mark gaps;
this notebook splits the belt at two of them, 2.5 AU (marked 3:1) and 2.82 AU (marked 5:2), into
the inner (blue), middle (orange) and outer (green) belt. Chart: Thunkii, after
Alan Chamberlain (JPL/Caltech),
Wikimedia Commons (CC BY 4.0).
# the main belt in three zones, split at the Kirkwood gaps at 2.5 and 2.82 AU
zone_edges = [2.1, 2.5, 2.82, 3.3]
zone_names = ["inner, 2.1-2.5 AU", "middle, 2.5-2.82 AU", "outer, 2.82-3.3 AU"]
in_belt = labelled[(labelled["a"] > 2.1) & (labelled["a"] < 3.3)].copy()
in_belt["zone"] = pd.cut(in_belt["a"], bins=zone_edges, labels=zone_names)
# the share of each kind in each zone
share_by_zone = pd.crosstab(in_belt["zone"], in_belt["complex"], normalize="index")
print(share_by_zone.round(2))
# three bars per zone, one per kind, each labelled with its share
zone_sizes = in_belt["zone"].value_counts().reindex(zone_names)
plt.figure(figsize=(8, 4.5))
for k, kind in enumerate(["S", "C", "X"]):
bars = plt.bar(np.arange(3) + (k - 1) * 0.27, 100 * share_by_zone[kind], width=0.27,
color=KIND_COLOURS[kind], label=kind)
plt.bar_label(bars, fmt="%.0f%%")
plt.xticks(range(3), [f"{name}\n{size} asteroids" for name, size in zip(zone_names, zone_sizes)])
plt.xlabel("zone of the main belt, nearest the Sun on the left")
plt.ylabel("share of the zone's classified asteroids (%)")
plt.ylim(0, 75)
plt.title(f"Stony S falls and carbon-rich C rises outward ({len(in_belt)} asteroids)")
plt.legend(title="kind")
plt.show()
Stony asteroids crowd the inner belt and carbon-rich ones the outer. S makes up 65% of the inner zone and 32% of the outer; C rises from 21% to 44%. This is the pattern Gradie and Tedesco found in 1982 (step 6 at the top), seen again in your 795 classified asteroids. So the orbit holds some information about the kind.
A random forest can use albedo and orbit together. You met it in the basalt week.
Random forest: Ask many slightly different trees and take a vote.
Each tree learns its own yes/no questions about
the columns, such as "is the albedo below 0.1?" and "is a above 2.8?", and the forest names the
kind most of its trees name.
How a random forest decides, in one picture:
Many decision trees each give an answer, and the forest takes their vote. Each tree is grown on its own bootstrap sample. Figure: TseKiChun (CC BY-SA 4.0), via Wikimedia Commons.
✏️ Your turn 4¶
How well does a forest on albedo and orbit sort the asteroids it learned from?
Why: the forest can combine albedo with where the asteroid orbits, which the rule cannot. If the orbit adds information, the forest should beat the rule.
Look back at: RandomForestClassifier, from Where was this basalt made? (10_where_basalt.ipynb).
four_inputs: the list of the four columnsalbedo,a,eandi.forest: aRandomForestClassifier(for examplen_estimators=100, random_state=0), fitted on those columns of every row oflabelled, withcomplexas the answer.forest_fit_f1: the macro F1 of the forest's predictions for the same rows oflabelled. Print it besiderule_f1.
Use these names, because the self-check looks for them: four_inputs, forest,
forest_fit_f1.
# ← your answer here: write your code here
# the self-check reads: four_inputs, forest, forest_fit_f1
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert sorted(four_inputs) == sorted(["albedo", "a", "e", "i"]), "the four columns named in the prompt"
assert hasattr(forest, "predict") and forest.n_features_in_ == 4, "forest is fitted on the four columns"
again = f1_score(labelled["complex"], forest.predict(labelled[four_inputs]), average="macro")
assert abs(forest_fit_f1 - again) < 1e-6, \
"forest_fit_f1: the macro F1 of the forest's kinds for the classified asteroids it was fitted on"
assert forest_fit_f1 > rule_f1, "on the asteroids it learned from, the forest should beat the two-cut rule"
print(f"✓ Your turn 4 — the forest scores {forest_fit_f1:.3f} on the asteroids it learned from")
On the asteroids it learned from, the forest scores 1.00, against 0.62 for the rule. That is too good to be true.
Overfitting: A curve that memorises the data you gave it fails on the data you did not.
A forest does the same: its trees can grow until every asteroid they saw sits in a leaf of its own, so a score on those same asteroids says nothing about the next one. The leaderboard asks about asteroids the model has never seen, so the next section asks the same.
3. How well does the classifier do on asteroids it has never seen?¶
A model is judged on data it did not learn from.
Train/test split: Hide some data from yourself, then check.
But which asteroids to hide? Copy how the leaderboard hid its own. Its 368 asteroids are 30% of the classified asteroids with a full colour spectrum, chosen at random so that each kind kept its share (a stratified split). Your hidden part must be made the way the leaderboard's was; this course calls that rule "split by structure". A split that is not random (all the big asteroids hidden, or one zone of the belt) would test a different question, and its score would not predict yours.
✏️ Your turn 5¶
How well does the forest do on asteroids it has never seen?
Why: this is the number that predicts how your model will do on the leaderboard. The score on the asteroids it learned from does not.
Look back at: train_test_split with stratify, from Where was this basalt made?
(10_where_basalt.ipynb).
- Split
labelledintofit_part(70%) andhidden_part(30%) withtrain_test_split, stratified bycomplex, withrandom_state=88. Take a.copy()of each, because you will add columns to them. Print the size of each, and how many of each kind the hidden part holds. - Fit the same kind of forest on
fit_partonly. Name itforest_before. - Add its kind for each hidden asteroid as the column
forestofhidden_part. test_f1: the macro F1 of that column against the true kinds. Print it besideforest_fit_f1andbaseline_f1.
Use these names, because the self-check looks for them: fit_part, hidden_part,
forest_before, test_f1.
# ← your answer here: write your code here
# the self-check reads: fit_part, hidden_part, forest_before, test_f1
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(fit_part) + len(hidden_part) == len(labelled) and 0.25 < len(hidden_part) / len(labelled) < 0.35, \
"the hidden part is about 30% of the classified asteroids, and every one is in exactly one part"
assert not set(fit_part["number"]) & set(hidden_part["number"]), "no asteroid in both parts"
again = f1_score(hidden_part["complex"], forest_before.predict(hidden_part[four_inputs]), average="macro")
assert abs(test_f1 - again) < 1e-6 and "forest" in hidden_part, \
"test_f1: the macro F1 on the hidden part of a forest fitted on the fit part, kept in the column forest"
assert baseline_f1 < test_f1 < forest_fit_f1 - 0.1, \
"the hidden score is out of range: was forest_before fitted on the fit part only?"
print(f"✓ Your turn 5 — the forest scores {test_f1:.3f} on the hidden asteroids, {forest_fit_f1:.3f} on those it learned from")
On asteroids it has never seen the forest scores 0.62, not 1.00. That is normal, and it is why the score on fitted asteroids is never the one to report. On the same hidden asteroids the two-cut rule, with the solution's cuts (0.075 and 0.16), scores about 0.59, so the forest is ahead. But the hidden part holds only 259 asteroids. A different 259 would give a slightly different score. How different?
Bootstrap: Ask the data the same question a thousand times, using a different random slice of itself each time.
Here a slice is a redraw of the hidden asteroids: as many as there are, picked at random with replacement, so one asteroid may come twice and another not at all. Each comes with its true kind and the forest's kind. A simple case first, on the six made-up asteroids:
# the six made-up asteroids, as a table: true kind and guessed kind
toy = pd.DataFrame({"complex": toy_truth, "guess": toy_guess})
# one redraw: six rows picked with replacement
toy_redraw = toy.sample(len(toy), replace=True, random_state=np.random.default_rng(2))
print(toy_redraw)
print("macro F1 of the six:", round(f1_score(toy["complex"], toy["guess"], average="macro"), 3))
print("macro F1 of this redraw:", round(f1_score(toy_redraw["complex"], toy_redraw["guess"], average="macro"), 3))
This redraw picked rows 5, 1, 0, 1, 2, 4: row 1 twice and row 3 not at all. Its macro F1 is 0.222, against 0.357 for the six as they are. Each row kept its own true kind and guess. A thousand redraws give a thousand scores, and their middle 95% is the interval:
Confidence interval: Not one number but the range your number would have wandered over, had the world rolled differently.
✏️ Your turn 6¶
How much could the test score move with other asteroids like these?
Why: two classifiers whose scores differ by less than this are not clearly different. You will need the interval to judge every comparison from here on.
Look back at: the bootstrap, from How long was the Earth frozen? (08_how_long_frozen.ipynb),
and the six made-up asteroids just above.
- Make a random generator,
rng = np.random.default_rng(88), and an empty list,redrawn_f1. - 1,000 times: redraw as many rows of
hidden_partas it holds, with replacement (random_state=rng), and append the redraw's macro F1 of itsforestcolumn against itscomplexcolumn. test_interval: the 2.5th and 97.5th percentiles of the 1,000 scores. Print it besidetest_f1.- Draw a histogram of the 1,000 scores, with dashed vertical lines at the two ends of the
interval and a solid one at
test_f1. A good figure puts the numbers on it: the legend says the interval and the test score, so a reader needs no printout.
Use these names, because the self-check looks for them: redrawn_f1, test_interval.
# ← your answer here: write your code here
# the self-check reads: redrawn_f1, test_interval
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert len(redrawn_f1) == 1000, "1,000 redrawn scores"
assert abs(np.mean(redrawn_f1) - test_f1) < 0.03, \
"each redraw's macro F1: its forest column against its complex column, over a redraw of the hidden part"
assert test_interval[0] < test_f1 < test_interval[1] and 0.03 < test_interval[1] - test_interval[0] < 0.25, \
"the 2.5th and 97.5th percentiles of the 1,000 redrawn scores, around the test score"
print(f"✓ Your turn 6 — 95% interval {test_interval[0]:.3f} to {test_interval[1]:.3f}")
The forest's test score is 0.62, with a 95% interval from 0.56 to 0.68. The interval is 0.12 wide. So a new classifier that scores, say, 0.02 more on these hidden asteroids has not clearly beaten the forest: other asteroids could reverse the order. One that scores above 0.68 would have.
One number hides which kind is missed. Your macro_f1 averaged three F1 scores; the next Your
turn looks at them one by one.
✏️ Your turn 7¶
Which kind does the forest miss?
Why: an average can hide one kind done badly. If one kind's F1 is far below the others, that kind is where a better classifier has to gain.
Look back at: your macro_f1 of Your turn 2.
f1_by_kind: a dict with keys"S","C"and"X", each holding that kind's F1 score on the hidden asteroids (theforestcolumn againstcomplex). Print it. The three average to yourtest_f1.- Draw it as a bar chart: one bar per kind in its colour, with its F1 written on the bar
(
plt.bar_label), a dotted line attest_f1, the y-axis from 0 to 1, axis labels and a title. The weak kind should stand out at a glance.
Use these names, because the self-check looks for them: f1_by_kind.
# ← your answer here: write your code here
# the self-check reads: f1_by_kind
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
by_kind = pd.Series(f1_by_kind, dtype=float)
assert sorted(by_kind.index) == ["C", "S", "X"], "f1_by_kind: one F1 score for each of S, C and X"
assert abs(np.mean(by_kind) - test_f1) < 1e-6, "the three F1 scores on the hidden asteroids average to test_f1"
assert by_kind.idxmin() == "X", "each kind's F1: that kind against the other two, on the hidden part's forest column"
print(f"✓ Your turn 7 — F1 on the hidden asteroids: S {by_kind['S']:.2f}, C {by_kind['C']:.2f}, X {by_kind['X']:.2f}")
S is sorted well (0.84), C fairly (0.71), and X badly (0.33). Of the hidden X asteroids, the forest called 26% X, 40% C and 34% S. Albedo and orbit place an X between the other two kinds, where both of them also live. To do better, the classifier needs something that tells X apart.
4. Does a different model, or more information, sort them better?¶
There are two ways to try to beat the forest. One is a different model that draws the boundaries between the kinds in another way. The other is more information: something about each asteroid that albedo and orbit do not hold. This section tries both, on the same hidden asteroids.
The SVM draws boundaries in a different way from a forest.
SVM: Of all the lines separating the two groups, take the one leaving the widest gap.
A forest cuts each column at a value, so its boundaries are steps; an SVM (with its default kernel)
draws smooth curved ones. Because it measures distances across the columns, an albedo of 0.1 and
an orbit of 2.8 AU must first be put on the same footing: StandardScaler rescales each column
to an average of 0 and a spread of 1, and make_pipeline(StandardScaler(), SVC()) joins the two,
as in the basalt week.
Of all the lines that separate two groups, the SVM takes the one with the widest gap; the points on the edge of the gap, the support vectors, decide where it lies. Figure: Larhmam (CC BY-SA 4.0), via Wikimedia Commons.
✏️ Your turn 8¶
Does an SVM on the same four columns beat the forest on the hidden asteroids?
Why: if the forest's step-shaped boundaries are what holds it back, a model with smooth boundaries should do better with the same information.
Look back at: make_pipeline, StandardScaler and SVC, from Where was this basalt made?
(10_where_basalt.ipynb).
svm: a pipeline that scales the columns and then fits an SVM, fitted onfit_part's four inputs andcomplex(the default settings are fine; others are too).- Add its kind for each hidden asteroid as the column
svmofhidden_part. svm_test_f1: its macro F1 on the hidden asteroids. Print it besidetest_f1andtest_interval; it should sit far abovebaseline_f1.
Use these names, because the self-check looks for them: svm, svm_test_f1.
# ← your answer here: write your code here
# the self-check reads: svm, svm_test_f1
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert hasattr(svm, "predict") and "svm" in hidden_part, \
"svm is a fitted model, and hidden_part has a column svm with its kind for each hidden asteroid"
again = f1_score(hidden_part["complex"], svm.predict(hidden_part[four_inputs]), average="macro")
assert abs(svm_test_f1 - again) < 1e-6, "svm_test_f1: the macro F1 of the svm column on the hidden part"
assert baseline_f1 + 0.2 < svm_test_f1 < 0.85, \
"the SVM's hidden score is out of range: was it fitted on the fit part's four columns, scaled?"
print(f"✓ Your turn 8 — the SVM scores {svm_test_f1:.3f} on the hidden asteroids, the forest {test_f1:.3f}")
The SVM scores 0.58, below the forest's 0.62 and inside its interval: a different shape of boundary bought nothing. Its F1 for X is 0.18. Both models are given the same two facts about an asteroid, how bright it is and where it orbits, and those facts do not separate X from its neighbours. So the next try is not a new model but more information.
Gaia measured each asteroid's colour, its reflectance spectrum (step 3 at the top). Does the spectrum tell X apart? The next cell draws the middle (median) spectrum of each kind.
# the middle (median) reflectance of each kind at each of the 14 wavelengths
wavelengths = [418, 462, 506, 550, 594, 638, 682, 726, 770, 814, 858, 902, 946, 990]
middle_spectrum = labelled.groupby("complex")[bands].median()
kind_sizes = labelled["complex"].value_counts()
plt.figure(figsize=(8.5, 5))
for kind in ["S", "C", "X"]:
plt.plot(wavelengths, middle_spectrum.loc[kind], "o-", color=KIND_COLOURS[kind], lw=2,
label=f"{kind}: {kind_sizes[kind]} asteroids")
# the three wavelengths the two colour numbers below are read at
for wavelength in [462, 770, 902]:
plt.axvline(wavelength, color="0.6", ls=":")
plt.text(wavelength, 0.73, f" {wavelength} nm", color="0.4")
# where olivine and pyroxene absorb, and the two numbers read from the S spectrum
plt.axvspan(800, 1000, color="tab:red", alpha=0.08)
plt.text(805, 1.25, "S dips here:\nolivine and pyroxene\nabsorb near 1,000 nm", color="tab:red")
plt.annotate("", xy=(770, middle_spectrum.loc["S", "R770"]), xytext=(462, middle_spectrum.loc["S", "R462"]),
arrowprops={"arrowstyle": "->", "color": "tab:red", "lw": 1.5})
plt.text(530, 1.17, "red slope:\nR770 minus R462", color="tab:red")
plt.xlabel("wavelength (nm): violet on the left, near infrared on the right")
plt.ylabel("reflectance, divided by its value at 550 nm")
plt.ylim(0.7, 1.35)
plt.title("The median spectrum of each kind: S is steep and dips, C is flat, X is in between")
plt.legend(loc="upper left")
plt.show()
print(middle_spectrum[["R462", "R770", "R902"]].round(3))
The three kinds have three shapes, and two numbers describe them.
- How red it is. S rises steeply from violet to red: at 770 nm it reflects
1.21 times its green light and at 462 nm only 0.84. C is nearly
flat, and X rises gently, between the two. The red slope,
R770minusR462, is 0.37 for S, 0.14 for X and 0.06 for C. - Whether it starts to dip. Only S falls after 800 nm, on its way down to the absorption
of olivine and pyroxene near 1,000 nm (step 3 at the top); Gaia's bands stop before the bottom
of that dip. The dip,
R770minusR902, is above zero only for S (0.045, against -0.048 for X and -0.034 for C).
These are the two things a planetary geologist reads first from such a spectrum. Neither is in albedo or orbit.
✏️ Your turn 9¶
Does the forest sort the hidden asteroids better when it also knows their colour?
Why: this is new information, not a new model. If colour records what an asteroid is made of, adding it should lift the score far beyond the interval, and above all for X.
Look back at: your Your turns 5 and 7, and the spectrum just above.
- Add two columns to both
fit_partandhidden_part:red_slope,R770minusR462, andband_depth,R770minusR902. six_inputs: the list offour_inputsand those two columns.colour_forest: the same kind of forest as in Your turn 5, fitted onfit_part's six inputs. Add its kind for each hidden asteroid as the columncolour_forestofhidden_part.colour_test_f1: its macro F1 on the hidden asteroids, andcolour_f1_by_kind: each kind's F1, as in Your turn 7. Print both besidetest_f1,test_intervalandf1_by_kind.
Use these names, because the self-check looks for them: six_inputs, colour_forest,
colour_test_f1, colour_f1_by_kind.
# ← your answer here: write your code here
# the self-check reads: six_inputs, colour_forest, colour_test_f1, colour_f1_by_kind
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
assert sorted(six_inputs) == sorted(["albedo", "a", "e", "i", "red_slope", "band_depth"]), \
"the four columns of Your turn 4 and the two colour columns"
assert np.allclose(hidden_part["band_depth"], hidden_part["R770"] - hidden_part["R902"]) \
and np.allclose(fit_part["red_slope"], fit_part["R770"] - fit_part["R462"]), \
"red_slope: the reflectance at 770 nm minus that at 462 nm; band_depth: at 770 nm minus at 902 nm"
again = f1_score(hidden_part["complex"], colour_forest.predict(hidden_part[six_inputs]), average="macro")
assert abs(colour_test_f1 - again) < 1e-6, \
"colour_test_f1: the macro F1 on the hidden part of the forest fitted on the fit part's six columns"
assert colour_test_f1 > test_interval[1], "with colour the forest should score above the interval of Your turn 6"
assert colour_f1_by_kind["X"] > f1_by_kind["X"], "each kind's F1 with colour, as in Your turn 7"
print(f"✓ Your turn 9 — with colour {colour_test_f1:.3f} (X {colour_f1_by_kind['X']:.2f}); without {test_f1:.3f} (X {f1_by_kind['X']:.2f})")
With two colour columns the forest scores 0.88 on the hidden asteroids, against 0.62 without: far above the top of the interval, 0.68. That is a clear win, from information rather than from a new model. The gain is largest for X, from 0.33 to 0.79: the gentle red slope with no dip is what marks an X.
The four classifiers on the same hidden asteroids. The next cell is given: it draws each one's macro F1 beside the forest's interval.
# hidden_part: Your turns 5, 8 and 9; rule_guess's cuts: Your turn 3
hidden_part["rule"] = np.where(hidden_part["albedo"] < low_cut, "C",
np.where(hidden_part["albedo"] < high_cut, "X", "S"))
classifiers = ["rule", "forest", "svm", "colour_forest"]
scores = []
for name in classifiers:
scores.append(macro_f1(hidden_part["complex"], hidden_part[name]))
print(pd.DataFrame({"classifier": classifiers, "macro_F1": scores}).round(3))
# what each classifier was given, for the labels on the figure
what_it_uses = ["two albedo cuts", "forest: albedo + orbit", "SVM: albedo + orbit",
"forest: albedo + orbit + colour"]
plt.figure(figsize=(8.5, 3.6))
plt.axvspan(test_interval[0], test_interval[1], color="0.85", label="the forest's 95% interval")
plt.axvline(baseline_f1, color="k", ls=":", label=f"always the most common kind: {baseline_f1:.2f}")
for k in range(len(classifiers)):
colour = "tab:green" if classifiers[k] == "colour_forest" else "0.35"
plt.plot(scores[k], k, "o", color=colour, markersize=10)
plt.text(scores[k] + 0.015, k, f"{scores[k]:.2f}", va="center", color=colour, fontweight="bold")
plt.yticks(range(len(classifiers)), what_it_uses)
plt.xlim(0, 1)
plt.xlabel("macro F1 on the hidden asteroids")
plt.ylabel("classifier: what it was given")
plt.title(f"Only colour lifts the score clear of the interval ({len(hidden_part)} hidden asteroids)")
plt.legend(loc="lower left")
plt.show()
Only colour clearly wins. The two-cut rule, the forest and the SVM, all built on albedo and orbit, sit within a few hundredths of each other; the forest with two colour columns sits far to the right of all three. The better classifier came from better information about what the surface is made of, not from a different algorithm.
One caution. The hidden asteroids have now been used to choose between classifiers, so the winner's score is a little kind to it: asteroids held back to choose between models are validation asteroids, and once they have chosen, they are no longer a fair test. The leaderboard, which nobody has used to choose anything, is the fair test of your final classifier.
The leaderboard¶
leaderboard holds 368 classified asteroids you have never seen. It has the input
columns of labelled (orbit, size, albedo and the colour bands) and row_id, the key of each
answer, but no name and no class. template has the same row_ids and an empty complex column
for your kinds.
- Your entry is your own classifier, fitted on the training file here; never look up a class online. Do not look up or download the hidden answers (the published classes of the leaderboard asteroids); an entry built from them is not your model, and it scores 0.
- The board shows a public score, the macro F1 on a random third of the 368 asteroids. The final score, on the other two thirds, is shown after the board closes, and your rank on it is graded. The board also shows a second final score on the asteroids smaller than 31.7 km.
- So choose your entry by your own fair test, not by the public score. The public third is about 123 asteroids; a classifier picked because it moved that number is tuned to them, and the rank uses the others.
- Up to 10 uploads that change
submission.csvare scored. An upload that changes only the notebook is free.
✏️ Your turn 10¶
Your leaderboard entry: what kind is each of 368 hidden asteroids?
Why: this is the project's answer, scored on classes you have not seen, and ranked against the other students. Choose by your fair test: the public score is a third of the asteroids, and your rank is decided on the rest.
Look back at: your Your turns 5 and 9.
- Take the classifier your fair test chose: the forest with the two colour columns, unless your
own tests on
fit_partandhidden_partfound something clearly better (all fourteen colour bands, an SVM on the colours, ...). Make a freshfinal_model, and name its list of input columnsfinal_inputs. - Add any made columns (
red_slope,band_depth) tolabelledand toleaderboard, made the same way. - Fit
final_modelon every row oflabelled: the hidden part is data too. submission: a copy oftemplate, with itscomplexcolumn filled with your classifier's kind for each leaderboard asteroid. Save it withsubmission.to_csv("submission.csv", index=False).- Print its first rows, its number of rows and of blanks, and how many of each kind it names.
- Run the self-check, then download
submission.csvand this notebook and upload both together (the steps are under "What you hand in" above) to Final project: asteroids.
Use these names, because the self-check looks for them: final_model, final_inputs,
submission.
# ← your answer here: write your code here
# the self-check reads: final_model, final_inputs, submission
# Self-check: run this after your answer. You do not need to read it — it only says whether something is off.
saved = pd.read_csv("submission.csv")
assert list(saved.columns) == ["row_id", "complex"], "submission.csv holds the template's two columns, row_id and complex"
assert (saved["row_id"].values == template["row_id"].values).all(), "one row per leaderboard asteroid, in the template's order"
assert saved["complex"].isin(["S", "C", "X"]).all(), "a kind for every leaderboard asteroid: S, C or X, in capitals"
assert (saved["complex"].values == submission["complex"].values).all(), "submission.csv holds the table submission"
assert saved["complex"].nunique() == 3, "the entry should name all three kinds"
assert final_model.n_features_in_ == len(final_inputs), "the final model is fitted on its own input columns"
assert len(final_model.predict(labelled[final_inputs])) == len(labelled), "the final model reads the columns of labelled"
assert f1_score(labelled["complex"], final_model.predict(labelled[final_inputs]), average="macro") > 0.8, \
"the final model, fitted on every classified asteroid, should sort them far better than the baseline"
print(f"✓ Your turn 10 — submission.csv: {len(saved)} kinds, {np.mean(saved['complex'] == 'X'):.0%} of them X")
The question, answered¶
Brightness and orbit sort rocky, carbon-rich and metal-bearing asteroids only part of the way: on asteroids it had never seen, a forest on albedo and orbit scored a macro F1 of 0.62, against 0.21 for always guessing S, and an SVM did no better. Two numbers read from each asteroid's colour, how red it is and whether it starts to dip after 800 nm, lifted the score to 0.88, most of all for the X group. What an asteroid is made of shows in the light it reflects, so the classifier improved when it was given that light.
The open question¶
Can a classifier trained on big asteroids be trusted on the small ones that are most of the belt, and can colour ever tell metal from rock?
Every asteroid you trained on was classified because it is big and bright enough for a telescope spectrum. Most asteroids are neither. The next cell, given, runs your forest with colour from Your turn 9 on every asteroid in the training file that Gaia gave a full colour spectrum, classified or not, and draws the share of each kind it names for the small unclassified ones beside the true shares of the classified ones.
# every asteroid in the training file with all 14 colour bands, classified or not
coloured = asteroids[asteroids[bands].notna().all(axis=1)].copy()
coloured["red_slope"] = coloured["R770"] - coloured["R462"]
coloured["band_depth"] = coloured["R770"] - coloured["R902"]
coloured["model_kind"] = colour_forest.predict(coloured[six_inputs]) # Your turn 9
# the classified ones, and the unclassified ones smaller than 10 km
small = coloured[coloured["complex"].isna() & (coloured["diameter"] < 10)]
print("asteroids with all 14 bands:", len(coloured), " classified:", coloured["complex"].notna().sum())
print("middle diameter, classified:", round(labelled["diameter"].median(), 1), "km",
" all with colour:", round(coloured["diameter"].median(), 1), "km")
print("kinds your model names for the", len(small), "unclassified asteroids under 10 km:")
print(small["model_kind"].value_counts(normalize=True).round(2))
# the share of each kind: true kinds of the classified asteroids, the model's kinds for the small unclassified ones
share_classified = labelled["complex"].value_counts(normalize=True)
share_small = small["model_kind"].value_counts(normalize=True)
plt.figure(figsize=(7.5, 4))
for k, kind in enumerate(["S", "C", "X"]):
left = plt.bar(k - 0.2, 100 * share_classified[kind], width=0.38, color=KIND_COLOURS[kind], edgecolor="k")
right = plt.bar(k + 0.2, 100 * share_small[kind], width=0.38, color=KIND_COLOURS[kind], edgecolor="k",
hatch="//", alpha=0.6)
plt.bar_label(left, fmt="%.0f%%")
plt.bar_label(right, fmt="%.0f%%")
plt.bar([0], [0], color="white", edgecolor="k", label=f"classified, true kind ({len(labelled)})")
plt.bar([0], [0], color="white", edgecolor="k", hatch="//", label=f"unclassified under 10 km, your model's kind ({len(small):,})")
plt.xticks(range(3), ["S", "C", "X"])
plt.xlabel("kind")
plt.ylabel("share of the asteroids (%)")
plt.ylim(0, 75)
plt.title("Do the small asteroids differ, or does the model stretch?")
plt.legend(loc="upper right")
plt.show()
The classified asteroids have a middle diameter of 31 km; the 31,769 asteroids with a full colour spectrum, 5 km. Only 175 of your 861 classified asteroids are smaller than 10 km, and they are the bright ones: their middle albedo is 0.27, against 0.21 for the unclassified ones that small. Your classifier still names a kind for every small one: for the 25,796 unclassified asteroids under 10 km it says 57% S and 20% X, against 45% and 23% among the classified ones. Is that what the small asteroids are, or the classifier stretched beyond the asteroids it learned from? Your fair test cannot say, because its asteroids were chosen the same way. After the board closes, the leaderboard's second final score, on its asteroids smaller than 31.7 km, gives a first hint. And whether an X like Psyche is metal, or rock that only has metal's colour, no colour can settle: that waits for the spacecraft.