In this homework, we will classify rock samples from a regional whole-rock geochemistry survey. The survey contains 10,000 samples. For each sample, a laboratory measured the concentrations of the major-element oxides (in weight percent), the bulk density, and the magnetic susceptibility. Field geologists mapped each sample site and assigned a lithology label: granite, basalt, or andesite. The classes are imbalanced.
The table is generated by the course package mlgeo_synth with a fixed seed; the instructor holds a hidden-seed variant used to spot-check submitted results.
In this homework we will train several classifiers to predict the class of a rock sample based on the measurements (features). We will practice data prep, dimensionality reduction, model design and training, model comparison, and feature importance selection.
Importing Libraries¶
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
1. Data Preparation (20 points)¶
We follow the following steps:
- read (1 point)
- clean (3 points)
- correlate (4 points)
- explore, spread of values (3 points)
- dimensionality reduction (9 points)
import mlgeo_synth
geo = mlgeo_synth.geochem_table(n=10000, seed=2026)
geo.insert(0, "sample_id", [f"S{i:05d}" for i in range(len(geo))])
geo.to_csv("rock_survey.csv", index=False)1.1 Data read¶
Read the pandas data frame from the csv file “rock_survey.csv”.
Task: read pandas data frame (1 point)
Save a copy of the data frame just in case.
Description of the data fields
- sample_id = Sample identifier, self explanatory.
The major-element oxides, reported in weight percent (wt%) of the bulk rock:
- SIO2 = silica (SiO2), the main oxide in most crustal rocks; high in felsic rocks like granite, low in mafic rocks like basalt.
- AL2O3 = alumina (Al2O3), hosted mainly in feldspars.
- FEO = total iron reported as FeO; high in mafic rocks.
- MGO = magnesia (MgO), hosted in olivine and pyroxene; high in mafic rocks.
- CAO = lime (CaO), hosted in calcic plagioclase and pyroxene.
- NA2O = soda (Na2O), hosted in sodic plagioclase.
- K2O = potash (K2O), hosted in alkali feldspar and mica; high in granite.
The oxides are subject to compositional closure: they sum to roughly 100 wt%, so when one goes up the others must go down.
The physical properties:
density_g_cm3 = bulk density of the sample in g/cm3.
mag_susc_si = magnetic susceptibility, in SI volume units; sensitive to the magnetite content of the rock.
label = lithology mapped by the field geologists (granite, basalt, or andesite). This will be the response variable which we will be trying to predict.
1.2 Data Cleaning¶
Basic stats about our dataset.
Task: Provide basic info for the pandas dataframe head (0.5 point)
Task: Find the data types of the database (floats, string, etc etc) using the info() function (0.5 point).
Are there any obvious feature (or element of the dataframe) that should not impact our prediction?
sample_id is just an identifier for accessing the rows back when they were stored in the original survey database. Therefore we will not need it for classification as it is not related to the outcome.
Task: Drop this column in the pandas dataframe. (1 point)
Find our how many examples there are, how many attributes or feature, and the type of class.
Task: How many objects are in each class? (1 point)
The classes are “granite”, “basalt”, and “andesite”. They are defined as strings, but we will convert them to integer in order to apply a loss function on the class labels during training. For this, we use the sklearn.preprocessing.LabelEncoder() function. We will do so and modify the classes in the dataframe. We should keep a copy of the original data frame to be safe.
1.3 Data correlations¶
Now let’s find the most basic correlations among features. This can be done using the corr() function to apply on the pandas dataframe. Evaluate this function and comment on what feature is correlated among others. It is convenient to use the matplotlib function matshow() for clarity. seaborn is a python module that makes really pretty statistical plots https://
Task: Plot the correlation matrix that can be called in the pandas dataframe. (2 points)
Hints:
Use functions of heatmap, add the labels in the axes. The colormap coolwarm is nice for divergent scales like correlations that vary between -1 and 1. The argument center=0 ensures that the colormap is divergent from zero. Make sure to ignore the label column “label”. Remember that dropping a column can be done in place rock_df.drop('label', axis=1).
Task: Reproduce the same plot for each of the three classes. (1 point) You can select the values from the pandas dataframe by selecting over the column ‘label’.
Task: Can you comment on groups of features that are correlated with each other or that appear independent from each other given these correlations? (1 point) You should expect strong correlations among the oxides: compositional closure forces them to trade off against each other, and magmatic differentiation drives them together. Do density and magnetic susceptibility correlate with the oxides? Do the patterns differ between the three lithologies?
1.5 Data exploration¶
Given the structure of the correlations, we will explore the values of the data.
1.5.a. Distributions of SiO2¶
Silica content is the first number a petrologist looks at: it increases with magmatic differentiation and separates felsic from mafic rocks.
Task: plot histograms for the ‘SIO2’ feature column for each class (1 point).
Task : Describe briefly the difference between the three histograms. (0.5 point)
Granite:
Basalt:
Andesite:
1.5.b. Density and magnetic susceptibility¶
We will now plot the bulk density (density_g_cm3) versus the magnetic susceptibility (mag_susc_si) colored by class. You can use the scatterplot or lmplot function in seaborn (https://
Task: do you see any obvious differences such that one could easily discriminate between the classes? (0.5 point)
1.5.c The major oxides¶
Recall: the correlation matrix shows that the major-element oxides are correlated with each other for all three classes.
Task: Plot histograms of the other oxides (AL2O3, FEO, MGO, CAO) and discuss why you expect these features to be correlated (1 points)
1.6 Data Dimensionality Reduction¶
At this point, we are left with 9 features: the seven oxides (SIO2, AL2O3, FEO, MGO, CAO, NA2O, K2O), density_g_cm3, and mag_susc_si. Among these, the oxides are correlated to each other. There is therefore a potential for reducing the dimensions of the features using PCA on these 7 features.
We will use the sklearn function sklearn.decomposition.PCA() to fit and transform the data into the PC coordinates. Let’s first explore how many PCs we need. Fit the PCA function over the total number of oxides. You will fit the PCA function over an array with the columns selected from the dataframe.
Task: Perform the PCA over a max number of PCs, output the explained variance ratio values, decide on an appropriate maximum number of PC to use (6 points)
Answer on how many PCs to use
We will now re-perform PCA with the number of PCs you found is most appropriate. Re-apply the fit-transform function. Update the dataframe by adding the PCA value(s) and dropping the columns of the 7 oxide features.
Task: PCA again, fit and transform, update the dataframe with the new feature(s) (3 points)
2. Unsupervised Clustering with KMeans (20 points)¶
In this section, we will explore if the data features will be sufficient for classification. As a first exploration, we will perform unsupervised classification with Kmeans clustering.
2.1 Perform preliminary Kmeans (10 points)¶
Implement Kmeans here for a given number of clusters and on the features of interest. Choose 3 features (for example PC1, density_g_cm3, and mag_susc_si; remember to scale them).
- Use
sklearnto perform Kmeans. - Repeat Kmeans and discuss (in a markdown cell) the stability of clustering (e.g., use visualization to qualitatively assess the stability).
2.2 Find the optimal number of clusters (5 points)¶
Use a method to establish the optimal number of clusters.
2.3 Discuss performance of clustering (5 points)¶
- Perform silhouette analysis (silhouette visualization and score)
- Calculate (python cell) and discuss (net markdown cell) homogeneity with respect to the ground truth labels using 3 appropriate metrics.
Question: After performing KMeans clustering and calculating the completeness, homogeneity, and Fowlkes-Mallows scores, how can you determine if these scores are good? Compare the obtained scores to the ideal values and explain what each score indicates about the clustering quality. What do you find from your results?
3 Machine Learning Models (30 points)¶
We will now train different models on this dataset. We have the features that remain after dimensionality reduction, 3 classes, and 10,000 samples. We will use K-Nearest Neighbors, Naive Bayes, Random Forest, Support Vector Machine, Histogram Gradient Boosting.
We now follow a normal machine learning workflow:
- Feature scaling (3) points)
- Train/test set split (2 points)
- Model design, training, testing (15 points)
- Model comparisons, pick your winner, discuss feature importance using Random Forest. (10 points)
3.1 Feature Scaling¶
Scaling all values to be within the (0, 1) interval will reduce the distortion due to exceptionally high values and make some algorithms converge faster. You can scale the features only by dropping the “label” column without modifying the dataframe in place, using the pandas function drop().
Task: Scale just the features (3 points)
3.2 Test, train, validation data sets.¶
Task: Split the data into a training and a test part. (2 points)
The models will be trained on the training data set and tested on the test data set. Use a stratified split (stratify=y) — the classes are imbalanced.
Computation time is important to account for when scaling up the data set and the model size. You can evaluate the relative computational time using the function time.perf_counter() to evaluate the absolute time. Then compare the computational time by making the difference between two time stamps:
t1=time.perf_counter()
t2=time.perf_counter()
tcomp = t2 - t1
We will also assess the model performance of these multi-class classifiers. We will evaluate the average of the scores over the 3 class labels.
In the following, we will be testing over several classifiers. Follow the steps:
- model definition/design
- training
- prediction on test
- evaluation: a) print the classification_report; b) save the precision, recall, fscore and accuracy in variables
3.3.a K Nearest Neighbors (3 points)¶
Check out the function arguments and definition here: https://
3.3.b Naive Bayes (3 points)¶
Check out the sklearn tutorial pages here: https://
Naive Bayes assumes the data to be normally distributed which can be achieved by scaling using the MaxAbsScaler. For this example then we will use the unscaled data, then rescale it.
3.3.c Random Forest Classifier (3 points)¶
Check out the tutorial page here: https://
3.3.d Support Vector Machine Classifier (3 points)¶
Check out the sklearn information page here: https://
3.3.e Histogram Gradient Boosting (3 points)¶
Check out the information page here: https://
3.4 Model performance and comparison¶
3.4.a Confusion Matrix and interpretation¶
Task: Plot the confusion matrix (2 points)
Use sklearn ``ConfusionMatrixDisplay" from the sklearn to visualize the confusion matrix
Task: Comment on what you see the best classifier is likely to be (1 point). You can also comment on the misclassification and confusion rates.
3.4.b K Fold Cross Validation¶
We will now perform k fold cross valdiation for the classifiers. We use the function cross_val_score on each estimator, on the training set, with 10 stratified folds (StratifiedKFold(n_splits=10)), and use macro-F1 as the score metric (scoring="f1_macro") — the classes are imbalanced, so accuracy would reward ignoring the rare class.
Task: perform the cross validation over K folds, output the mean and standard deviation of the macro-F1 score (3 points)
Task: Which method won the Xval test (1 point) ?
see the cell below
3.4.c And the winner is ...¶
Let’s compare the results. Task: Create a pandas dataframe with all of the performance metrics, including the results from K-fold cross validation. (2 points)
Task: Comment on the macro-F1 scores and performance and choose a winner. (1 point)
see the cell below
4 Summary (4 points)¶
4.1 Feature Importance using Random Forest Classifier¶
Decision Trees have the unique property of being able to order features by their ability to split between the classes. If some features dominate over other in the predictive power of classes, one can further reduce the dimension of the features for additional analysis. The vector of feature importance is the module rfc.feature_importances_, sorted with ascending importance. Store the vector of importance.
Recall the caveat from lesson 3.7: importance is not causation. The ranking reports what the model uses to predict, not what makes a rock a granite — and impurity importances split credit among correlated features, including the PCs you built from the oxides.
Task: plot a bar plot using the function matplotlib.pyplot.bar. (2 points)
Task: What are the top three features (1 point)?
enter in the cell below
In this notebook, you have likely found that the differentiation-related features (the first oxide PC, density, magnetic susceptibility) separate the lithologies. A petrologist would have told you that silica content separates granite from basalt — but now you can quantify how well, and what it costs to automate.
Task: Briefly comment on what you have learned (1 point)
see the cell below