Gaussian Anomaly Detection

Anomaly detection is one of the most practical applications of unsupervised and semi-supervised machine learning. Please note dear reader that the next chapter introduces K-means that is another useful semi-supervised machine learning.

Rather than asking “what class does this belong to?”, we ask a subtler and often more useful question: “does this observation look like the ones I have already seen?” In fraud detection, manufacturing quality control, medical diagnostics, network intrusion monitoring, and predictive maintenance, the interesting events are usually rare. Trying to build a supervised classifier for events that occur one time in a thousand often fails: there are simply not enough positive examples to learn a good decision boundary. Anomaly detection sidesteps this problem entirely by modeling the normal distribution of data and flagging observations that fall outside it.

In this chapter we implement a multivariate Gaussian anomaly detector in Gerbil Scheme from scratch, then apply it to the classic Wisconsin Breast Cancer dataset. The detector fits a Gaussian distribution to each feature independently (a “diagonal covariance” assumption), computes the probability density of each observation, and flags observations whose density falls below a threshold as anomalies. The threshold is selected automatically by sweeping over 200 candidate values and picking the one that maximizes performance on a cross-validation set.

The core idea is that if we assume each feature of a “normal” observation is drawn from a Gaussian distribution with mean and variance , then the probability density of a full observation vector can be approximated by combining the per-feature Gaussians. The Gaussian probability density function for a single feature is:

A standard multivariate Gaussian anomaly detector multiplies these per-feature densities together (assuming feature independence). The implementation in this chapter uses a slight variation: it computes the arithmetic mean of the per-feature PDFs rather than their product. This has the practical benefit of keeping the resulting numbers on a comfortable scale (avoiding vanishingly small products when there are many features) while still giving a monotone score suitable for thresholding.

The training procedure is straightforward:

1. Split the examples into training (_{60%), cross-validation (}20%), and test (~20%) sets. The training set is filtered to contain mostly normal examples.
2. Fit the per-feature mean and variance using only the training set.
3. Sweep many candidate values of and pick the one that yields the fewest misclassifications on the cross-validation set.
4. Evaluate the final model on the held-out test set, reporting precision, recall, and -score.

A key trick used here (following Andrew Ng’s classic Coursera Machine Learning lectures) is that we deliberately let a small number of anomalous examples leak into the CV and test sets - they are essential for tuning and for measuring recall.

We use the Wisconsin Breast Cancer Dataset, a well-known benchmark in medical machine learning. Each row is derived from a fine needle aspirate (FNA) of a breast mass and contains 9 cytological features scored on an integer scale from 1 to 10, plus a target label:

• Features (columns 1-9): clump thickness, uniformity of cell size, uniformity of cell shape, marginal adhesion, single epithelial cell size, bare nuclei, bland chromatin, normal nucleoli, mitoses.
• Target (column 10): 2 for benign, 4 for malignant.

The raw CSV format is refreshingly simple. Here are the first few lines of data/cleaned_wisconsin_cancer_data.csv:

The final column of 2 or 4 is what we want to remap to {0.0, 1.0} so the detector can distinguish anomalies (malignant) from normals (benign).

The project directory source_code/anomaly-detection contains:

This module contains the entire anomaly detection engine. It has no dependency on the Wisconsin dataset - it can be reused for any labelled dataset in which the final column encodes the anomaly label as 1.0 (anomaly) or 0.0 (normal).

We start with the imports and exports. The :std/format module gives us format for printing, :std/pregexp is available for regex use in downstream modules, and :gerbil/gambit exposes the underlying Gambit runtime primitives.

The constant SQRT_2_PI is precomputed once as because it appears in the denominator of every Gaussian PDF evaluation. Note the small trick (acos -1.0) to obtain without hardcoding a magic number.

Gerbil’s default number printing yields long decimal expansions (like 0.8571428571428571), which are inconvenient for reporting metrics. The next two helpers give us format-decimal for pretty-printing floats to a fixed number of places:

format-decimal multiplies by , rounds to the nearest integer, divides back, then post-processes the string to zero-pad the fractional part and prepend a leading 0 when necessary. The result is a stable, readable representation like "0.8571".

The split-data procedure randomly divides the examples into three sets. Crucially, it applies an additional filter to the training set: it accepts an example into training only if it is a normal one (target < 0.5), or with a small 10% probability if it is an anomaly. This mirrors the standard anomaly detection recipe - train primarily on normals, but keep some anomalies in CV and test for threshold tuning and evaluation.

The function returns an association list keyed by 'train, 'cv, and 'test. Each example is a vector where index (- num-features 1) holds the target label.

Fitting is textbook: compute the sample mean and sample variance for each feature. compute-mu loops over all examples, accumulating a per-feature sum into a vector, then divides by the number of examples:

compute-sigma-sq follows the same shape but accumulates squared deviations from the mean. Note the small floor of 1.0e-10 to prevent divide-by-zero when a feature happens to be perfectly constant on the training set:

Notice that the outer loop terminates at (- nf 1) rather than nf. This deliberately skips the last column (the target label) - we don’t want to include the label in the Gaussian model.

The gaussian-p function computes the density score used for classification. As mentioned above, it takes the arithmetic mean of the per-feature Gaussian PDFs rather than their product. This is the score we will threshold with :

build-detector is the constructor. It takes raw labelled examples, splits them, computes , and stashes everything into a hash table that carries the model state:

The core of training is to try many candidate values of and pick the one that produces the fewest classification errors on the cross-validation set. train-helper recomputes (from the training set) and counts CV errors for a given epsilon:

The subtle bit is the inner if: for a true anomaly, an error occurs when the score is above (we called it normal); for a true normal, an error occurs when the score is below (we called it anomalous).

test-model runs the full confusion matrix on the held-out test set and prints precision, recall, and :

Finally, train sweeps 200 values of starting from 0.001 in increments of 0.005, keeps the value that minimized CV errors, then evaluates on the test set:

is-anomaly is the user-facing predicate: given a trained detector and a new observation vector, return #t if the Gaussian score falls below the learned threshold.

The demo file glues the detector to the Wisconsin dataset. It imports the detector module (via the sibling filename import "detector") plus a few standard modules for I/O and regex splitting:

load-csv reads the file line-by-line, trims whitespace, splits on commas, converts each token to a number, and returns a list of vectors:

Empty lines are silently skipped, which makes the loader forgiving of trailing newlines.

Raw integer scores from 1 to 10 are not ideal for Gaussian modeling - the features are heavily skewed, with many observations pinned at 1. The preprocessing pipeline in preprocess-wisconsin applies three transformations to each row:

1. Scale down the raw features (indices 0-8) by 0.1, moving them into a smaller numeric range.
2. Log-transform each feature as . The additive constant prevents and softens the tail.
3. Per-row min-max scale the features to the interval , so no feature dominates the Gaussian by virtue of its raw magnitude.

The target column is separately remapped from to :

The main function ties everything together: load the CSV, preprocess, build the detector, train it (which internally sweeps and evaluates on test), then run a couple of light sanity checks:

Run either make test or invoke the interpreter directly. The Makefile is minimal:

A typical run produces output similar to the following. Because split-data uses random sampling, the exact numbers will vary from run to run:

The output reports several important quantities.

• Loaded 648 examples: the total number of records read from the CSV.
• Training / Cross-val / Test set sizes: notice that the training set is smaller than a naive 60% split would give. That is expected - split-data deliberately drops most of the anomalous rows from the training set so the model learns what “normal” looks like.
• Best epsilon = 0.8910: the classification threshold selected by the sweep. Any observation whose Gaussian score falls below 0.8910 is flagged as an anomaly.
• TP=24 FP=4 FN=1 TN=39: the confusion matrix on the 68-observation test set. TP (true positive) is a correctly-identified malignant case, FN (false negative) is a missed malignant case, FP (false positive) is a benign case wrongly flagged, and TN (true negative) is a correctly-identified benign case.
• precision = 0.8571: of the 28 observations the model called “anomaly”, 24 were truly malignant. Precision answers “when the model raises an alarm, how often is it right?”
• recall = 0.9600: of the 25 truly malignant observations, the model caught 24. Recall answers “of all the actual anomalies, how many did we catch?”
• F_1 = 0.9057: the harmonic mean of precision and recall. F_1 is a single-number summary of overall detection quality when both false positives and false negatives matter.

The recall value is particularly encouraging in a medical setting - missing a malignant tumor is a much more costly mistake than a false alarm, so a recall above 0.95 is a genuinely useful outcome.

We built a full-featured Gaussian anomaly detection engine in Gerbil Scheme in fewer than 200 lines of code. The key ideas transfer to any application where “the interesting event is rare”:

• Fit a Gaussian per feature on your normal data.
• Use the density function as a score.
• Use a held-out cross-validation set with a few real anomalies in it to pick the best classification threshold .
• Report precision, recall, and on a completely held-out test set.

Gerbil Scheme turns out to be a comfortable environment for this kind of numerical work. Vectors give us cheap random-access storage, hash tables give us a natural home for the trained model’s state, and the top-level def / import / export module system keeps detector.ss cleanly reusable across different datasets. Try pointing wisconsin_demo.ss at your own CSV data - as long as the last column is a binary anomaly label, the detector will happily consume it.

1. Product-of-PDFs Score: gaussian-p currently returns the arithmetic mean of the per-feature PDFs. Add a new function gaussian-p-product that instead returns the classical multivariate Gaussian score (the product of the per-feature PDFs) and compare its precision/recall/F1 against the mean-based score on the Wisconsin dataset.

2. Sensitivity to the Preprocessing Pipeline: In preprocess-wisconsin, remove the log-transform step (leave only the *0.1 scaling and the min-max normalization). Re-run the demo and report how the best , precision, recall, and change. Explain why the log-transform helps or hurts.

3. Configurable Sweep Range: Modify train in detector.ss to accept optional keyword arguments eps-start, eps-step, and eps-count so users can control the epsilon sweep without editing the source. Provide sensible defaults matching the current behavior (0.001, 0.005, 200).

4. Save and Load a Trained Detector: Implement save-detector and load-detector functions in detector.ss that serialize the trained , , best-eps, and num-features to a plain-text file and read them back in a later session, so that training does not need to be repeated to classify new observations.

5. Try a Different Dataset: The KDD Cup 1999 network intrusion dataset is a classic benchmark for anomaly detection. Write a kdd_demo.ss that mimics wisconsin_demo.ss for a subset of the KDD data, including any needed preprocessing (categorical features must be converted to numeric form). Report how well the detector performs.