library("bayesplot")
library("cmdstanr")
library("ggplot2")
library("khroma")
library("quartoExtra")
# Source for the modified reliability plot
source("../../code/helpers.R")
theme_set(good_theme)
bayesplot_theme_set(good_theme)
color_scheme_set(scheme = c(unname(colour("vibrant")(7)[c(3,7,4,1,2,5)])))
scale_colour_discrete = scale_colour_vibrant
scale_fill_discrete = scale_fill_vibrant
source("../../code/helpers.R")
SAVE_MODEL = TRUECalibration plots for the easy tasks of identifying penguin species in the Palmer Penguins data set.
For the set of \(N\) observations of \(K\) features \(x_{k,n}\) and their corresponding species \(y_{n} \ in \{1,\dots,D\}\), we perform the classification using the categorical logit classifier, that is using the following likelihood \[ p(y_n = s_j \mid \mathbf{w}) = \frac{\exp\left(\beta_{n,j}\right)}{\sum_{d=1}^D\exp\left(\beta_{n,j}\right)}, \] where \[ \beta_{n,j} = w_{0,j} + \sum_{k=1}^Kw_{k,j}x_{k,n} \]
data {
int N; // number of observations
int D; // number of features
int N_classes; // number of classes
matrix [N, D] X; // observation data
array[N] int <lower = 1, upper = N_classes> y; // target values {1,..., N_classes}
}
transformed data {
// Normalize the data
matrix[D + 1, N] X_stn;
X_stn[D + 1, ] = rep_row_vector(1, N);
for (d in 1:D) {
X_stn[d,] = to_row_vector((X[, d] - mean(X[, d])) / sd(X[, d]));
}
}
parameters {
// Matrix of the linear weights.
matrix[N_classes, D + 1] W;
}
transformed parameters {
// Compute the linear transformations
matrix[N_classes, N] logits;
for (c in 1:N_classes) {
logits[c, ] = W[c, ] * X_stn;
}
}
model {
// Evaluate prior density
for (d in 1:(D + 1)) {
for (c in 1:N_classes) {
target += normal_lpdf(W[c, d] | 0, 1);
}
}
// Evaluate likelihood
for (n in 1:N) {
target += categorical_logit_lpmf(y[n] | logits[,n]);
}
}
generated quantities {
vector[N] yrep; //posterior predictive samples
matrix[N,N_classes] lpd; //log posterior predictive densities
for (n in 1:N) {
yrep[n] = categorical_logit_rng(logits[,n]);
}
for (n in 1:N) {
for (c in 1:N_classes) {
lpd[n, c] = categorical_logit_lpmf(c | logits[,n]);
}
}
}After loading the model, we pass the data and run the inference using HMC. As a result, we obtain 4000 samples from the posterior as well as posterior predictive samples and posterior predictive densities.
Running MCMC with 4 parallel chains...
Chain 4 finished in 6.0 seconds.
Chain 3 finished in 6.1 seconds.
Chain 2 finished in 6.3 seconds.
Chain 1 finished in 6.5 seconds.
All 4 chains finished successfully.
Mean chain execution time: 6.2 seconds.
Total execution time: 6.7 seconds.The posterior distributions of the weights show reasonable separation between classes (the first index). The dimension with most overlap seems to be the intercept term, \(W[j,5]\) below.
The common approach of plotting a bar char of the observations overlaid with posterior means and 95% confidence intervals only gives a crude idea of the calibration of the model predictions.
From the following three reliability diagrams, one can see, that although sometimes underconfident in the predicted probabilities, the model is capable of very clearly separating the three species.
