spine

Partition data into scale strata. One call: raw data in, strata out.
You supply a measure and an orientation. Those are the scientific input, and neither has a default. Everything after them is computed.

Syntax and options

Read from the function itself, so this can never disagree with help spine.
help spine
SPINE Apply SPINE to raw data in one call. [scores, Result] = spine(data, measure, orientation) [scores, Result] = spine(data, measure, orientation, 'Name', Value, ...) [scores, Result, Info] = spine(...) MEASURE and ORIENTATION are both required. They are the scientific input: SPINE does not know what your data means, and choosing them is the modelling step. Neither has a default. Accepted data numeric vector one value per object numeric matrix rows are objects, columns are features struct array one element per object; numeric scalar fields become the features, a text field becomes the identifiers table numeric variables become the features, row names or a text variable become the identifiers Measure Any name accepted by spineMeasure. Type "help spineMeasure" for the menu. Orientation "upper" large measure is atypical "lower" small measure is atypical "bilateral" atypical on either side of an anchor object Options 'Objects' "rows" (default) or "columns", for numeric matrices 'K' neighbour count for the knn measure (default 5) 'ObjectIds' identifiers, overriding anything found in the data Any remaining option is passed to spineOptions, so BoundaryMode, ProfileWeights, and BilateralAnchor all work here unchanged. Outputs scores one score per object, in YOUR object order Result the readable record, everything in your object order (below) Info the reference record returned by spineCore, left untouched Result fields scores one score per object, 0 typical to 1 most atypical data the ingested data, rows as objects, for plotting mu the computed measure, one value per object stratum stratum index per object nStrata number of strata nBoundaries number of cuts, always nStrata - 1 boundaryValues measure value at each cut, largest value below the cut boundaryIntervals the open interval each cut lies in, one row per cut strataRanges measure range of each stratum, one row per stratum objectIds the identifiers used measure what spineMeasure did: name, parameters, dropped columns orientation the orientation you chose order, inverse the sort mapping applied before calling spine. When objects were omitted, order holds the caller's row numbers of the kept objects in ascending measure order, and inverse is NaN at the omitted rows. Info the reference record, for spinePlot and for auditing boundaryValues is the transferable result. Info.boundaries are ranks into an internally sorted copy and mean nothing outside this call; the measure VALUE at a cut is what you apply to new data or quote in a paper. Because the measure is sorted before spineCore sees it, the index-valued fields of Info (permutation, sortedIds, boundaryLeftObjectIds, boundaryRightObjectIds) refer to that sorted copy, not to your ordering. Result holds the remapped view. This file sits at the top of the MATLAB folder and needs no setup. Calling it puts core and interface on the path for you. It calls SPINE only through spineCore(mu, orientation, options) and spineOptions. No file under core is read or modified. Author Abd AlRahman AlMomani Assistant Professor of Data Science and Mathematics, ERAU, AZ almomaniar@gmail.com

A first call

Two hundred objects from a standard normal cloud, and twelve sitting a long way out. Nothing tells SPINE that there are two groups.
rng(11);
X = [randn(200,2); randn(12,2) + 9];
 
[scores, Result] = spine(X, 'dist2mean', "upper");
spineSummary(Result)
SPINE 212 objects, 2 strata measure dist2mean | orientation upper stratum objects measure range score 1 200 0.04007 to 3.793 0.000 2 12 10.91 to 13.82 1.000 cut at 3.793 Result.boundaryValues applies these cuts to new data.
scores runs from 0 for the typical stratum to 1 for the most atypical, one value per object, in your object order.
spineDataPlot(Result, 'Title', 'a bulk and a far group');

How it works, in one picture

SPINE sorts your measure, then walks up it accumulating an incremental nested entropy: at each rank, how much uncertainty the prefix of objects seen so far carries. Adding another object of the same scale increases it. Adding an object from a genuinely larger scale does not: it concentrates the distribution, and the profile turns over.
A cut is placed where the profile turns over. That is the whole idea. There is no threshold and no target number of groups; the entropy-supported transition either happens or it does not.
rng(11);
X = [randn(200,2); randn(12,2) + 9];
[~, Result] = spine(X, 'dist2mean', "upper");
 
I = Result.Info;
figure('Color','w','Position',[100 100 620 300]);
plot(I.E, 'LineWidth', 1.4, 'Color', [0 0.35 0.7]); hold on
plot(I.candidates, I.E(I.candidates), 'o', 'MarkerSize', 8, ...
'MarkerEdgeColor', [0.85 0.33 0.10], 'LineWidth', 1.4);
plot(I.boundaries, I.E(I.boundaries), 'o', 'MarkerSize', 8, ...
'MarkerFaceColor', [0.85 0.33 0.10], 'MarkerEdgeColor', [0.85 0.33 0.10]);
xlabel('rank in the sorted measure'); ylabel('incremental nested entropy');
title('the profile rises with the bulk, then turns over where the scale changes');
legend({'profile', 'candidate transition', 'retained cut'}, ...
'Location', 'southeast', 'Box', 'off');
grid on; hold off
The turnover sits at rank 200 of 212, which is exactly where the twelve planted far objects begin. Nothing told SPINE there were two groups, or where to look.
fprintf('cut placed at rank %d of %d\n', I.boundaries, numel(I.E));
cut placed at rank 200 of 212

The cut is the transferable result

boundaryValues is a value in your own measure's units, so it can be applied to new data or quoted in a paper. Info.boundaries, by contrast, are ranks into an internally sorted copy and mean nothing outside the call.
Result.boundaryValues
ans = 3.7928

Any container you already have

A numeric vector or matrix, a struct array, or a table. Identifiers carried by a table or struct come back in Result.objectIds.
rng(8);
sites = [randn(55,3); randn(5,3) + 7];
T = array2table(sites, 'VariableNames', {'x','y','z'}, ...
'RowNames', compose("site%02d", (1:60).'));
 
[sT, ResT] = spine(T, 'robustrownorm', "upper");
disp(ResT.objectIds(sT > 0).')
"site56" "site57" "site58" "site59" "site60"

Missing values

Real data has gaps. The default refuses rather than guessing.
Xm = X; Xm(4,2) = NaN; Xm(150,1) = NaN;
try
spine(Xm, 'dist2mean', "upper");
catch err
disp(err.message)
end
X must be finite (no NaN or Inf).
'Missing', "omit" drops those objects and reports them. Their scores come back as NaN, so the result still lines up with your rows: row 150 is still row 150.
[sm, Resm] = spine(Xm, 'dist2mean', "upper", 'Missing', "omit");
fprintf('omitted %s, scores returned %d for %d input rows\n', ...
mat2str(Resm.omitted.'), numel(sm), size(Xm,1));
omitted [4 150], scores returned 212 for 212 input rows

Very large collections

A profile that cannot be resolved at binary64 precision is refused rather than approximated. This follows the particular values, not the object count, so there is no size threshold to stay under. When it happens, 'Thin', true locates the cut on a deterministically thinned copy and applies it to every object.
[scores, Result] = spine(mu, 'asis', "upper", 'Thin', true);
Result.thinned records how many values were retained.

See also

spineSummary for what a partition found, spineCompare to check whether a different measure would have said something different, spineMeasure for the ten named measures, and spineDataPlot and spinePlot for the figures.
close all