API — ICoreBlocks/ICoreMath/Statistics
The public contract of 18 header(s) under src/ICoreBlocks/ICoreMath/Statistics — 18 class/struct definition(s), 24 declaration(s). Each section shows the header's banner and its public (and protected-virtual) surface exactly as the file writes it.
ICoreClassificationScoring.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreClassificationScoring.h
ICoreClassificationScoring#
ICoreClassificationScoring.h:41 · class · 2 declaration(s)
Toolbox row T3.21 -- classify, confusionmat and perfcurve: the one classifier the Statistics Toolbox ships WITHOUT an object, and the two functions that score any classifier's output.
class ICoreClassificationScoring {
public:
// MATLAB's five discriminant types.
enum class Discriminant { Linear, Quadratic, DiagLinear, DiagQuadratic, Mahalanobis };
static bool discriminantOf(const std::string& word, Discriminant& type);
// `[class, err, posterior, logp] = classify(sample, training, group[, type[, prior]])`.
// `prior` may be empty for MATLAB's default (equal over the non-empty
// groups); `posterior` and `logDensity` come back empty for the
// Mahalanobis type, which is what MATLAB answers.
//
// `group` is numeric here: MATLAB's grouping variable may also be a cell
// of names, and this console has no cell kind (Z3). The group LEVELS are
// the distinct values sorted ascending, which is `grp2idx`'s own order.
static bool discriminantClassify(const ICoreMatrix& sample, const ICoreMatrix& training,
const ICoreMatrix& group, const Discriminant& type,
const ICoreMatrix& prior, ICoreMatrix& classes,
double& resubstitutionError, ICoreMatrix& posterior,
ICoreMatrix& logDensity, std::string* whyNot = nullptr);
// `[C, order] = confusionmat(actual, predicted[, "Order", order])`. `order`
// in is empty for the default; `groups` out is the list the rows and
// columns are in.
static bool confusionMatrix(const ICoreMatrix& actual, const ICoreMatrix& predicted,
const ICoreMatrix& order, ICoreMatrix& counts,
ICoreMatrix& groups, std::string* whyNot = nullptr);
// The performance criteria `perfcurve` can put on either axis. Each is a
// ratio of the four counts at one threshold, and the names are MATLAB's
// own (with its synonyms: `sens`, `reca` and `tpr` are one criterion).
enum class Criterion { TruePositiveRate, FalsePositiveRate, FalseNegativeRate,
TrueNegativeRate, Precision, NegativePredictiveValue, Accuracy };
static bool criterionOf(const std::string& word, Criterion& criterion);
// `[X, Y, T, AUC] = perfcurve(labels, scores, posclass[, "XCrit", x, "YCrit", y])`.
//
// The area is the trapezoid rule over the curve as returned, with any
// segment whose end is undefined (a precision at zero predicted positives)
// left out -- which is how MATLAB answers a number for a curve that starts
// at NaN.
static bool performanceCurve(const ICoreMatrix& labels, const ICoreMatrix& scores,
const double& positiveClass, const Criterion& horizontal,
const Criterion& vertical, ICoreMatrix& x, ICoreMatrix& y,
ICoreMatrix& thresholds, double& area,
std::string* whyNot = nullptr);
};
};ICoreClusterAnalysis.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreClusterAnalysis.h
ICoreClusterAnalysis#
ICoreClusterAnalysis.h:55 · class · 3 declaration(s)
The clustering zone: the hierarchical block of toolbox row T3.20 -- linkage, cluster, clusterdata, cophenet, inconsistent and silhouette -- and T3.18's kmeans, at the bottom of the file.
class ICoreClusterAnalysis {
public:
// MATLAB's seven method words. `Ward` is spelled `ward` on the console the
// way MATLAB spells it, though `linkage.m`'s own table writes it `ward's`.
enum class Method { Single, Complete, Average, Weighted, Centroid, Median, Ward };
static bool methodOf(const std::string& word, Method& method);
// Whether this method's tree is re-sorted by `rearrange` -- true for
// average, weighted, complete and Ward's, false for the other three.
static bool methodIsRearranged(const Method& method);
// `linkage(Y[, method])` over a CONDENSED distance row, T3.19's `pdist`
// order. Answers the (m-1) x 3 tree.
static bool linkageOfDistances(const ICoreMatrix& condensed, const Method& method,
ICoreMatrix& tree, std::string* whyNot = nullptr);
// `linkage(X[, method[, metric]])` over a data matrix: `pdist` first, then
// the tree. `metric` is one of T3.19's twelve words.
static bool linkageOfPoints(const ICoreMatrix& points, const Method& method,
const std::string& metric, const double& exponent,
ICoreMatrix& tree, std::string* whyNot = nullptr);
// `inconsistent(Z[, depth])` -- the (m-1) x 4 table: the mean and the
// standard deviation of the link heights in the sub-tree `depth` levels
// below each link, how many links that was, and the inconsistency
// coefficient itself.
//
// The standard deviation divides by n-1 EXCEPT when n is 1, where it
// divides by 1 and answers zero -- `inconsistent.m`'s `s(3)-(s(3)~=1)`,
// which is not the same as MATLAB's `std` and would be a NaN if it were.
static bool inconsistentTable(const ICoreMatrix& tree, const size_t& depth,
ICoreMatrix& out, std::string* whyNot = nullptr);
// Which criterion a cut is measured against.
enum class Criterion { Inconsistent, Distance };
// `cluster(Z, "Cutoff", c[, "Criterion", crit][, "Depth", d])` -- every
// link at or below the cutoff stays connected, and a link counts as
// connected only if every link beneath it is too (`checkcut` in
// cluster.m). Answers one cluster number per leaf, renumbered by first
// appearance.
static bool clusterByCutoff(const ICoreMatrix& tree, const double& cutoff,
const Criterion& criterion, const size_t& depth,
ICoreMatrix& out, std::string* whyNot = nullptr);
// `cluster(Z, "MaxClust", k)` -- at most k clusters.
//
// ⚠ THE DEFAULT CRITERION CHANGES WITH THE FORM. `cluster(Z, 3)` and
// `cluster(Z, "MaxClust", 3)` cut by DISTANCE, not by inconsistency, even
// though "inconsistent" is the documented default: `cluster.m` overrides
// it whenever a cutoff was not given and the criterion was not named. The
// inconsistency-based count form is reachable only by naming the
// criterion, and is what `inconsistentCount` is.
static bool clusterByCount(const ICoreMatrix& tree, const size_t& maxClusters,
ICoreMatrix& out, std::string* whyNot = nullptr);
// `cluster(Z, "MaxClust", k, "Criterion", "inconsistent"[, "Depth", d])`:
// the largest inconsistency cutoff that still leaves at most k clusters,
// found by walking the distinct coefficients downward.
static bool clusterByInconsistentCount(const ICoreMatrix& tree, const size_t& maxClusters,
const size_t& depth, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `cophenet(Z, Y)` -- the cophenetic correlation coefficient, and the
// cophenetic distances themselves in `pdist` order.
static bool cophenetic(const ICoreMatrix& tree, const ICoreMatrix& condensed,
double& coefficient, ICoreMatrix& distances,
std::string* whyNot = nullptr);
// `silhouette(X, idx[, metric])` -- one value per row of X.
//
// ⚠ THE DEFAULT METRIC IS SQUARED EUCLIDEAN, not Euclidean. `silhouette.m`
// says `distType = 'sqeuclidean'` when no metric is given, which is the
// one place in the Statistics toolbox where a squared distance is the
// default, and it moves every number.
//
// A point alone in its cluster has no within-cluster distance to average.
// MATLAB answers its silhouette as exactly 1 -- measured against R2026a,
// `silhouette([0 0; 1 1; 10 10], [1;1;2])` ends in a 1 -- because the
// average it cannot take comes out 0 and `(b - 0)/max(0, b)` is 1. Not a
// NaN, and not the 0 a "no information" reading would give.
static bool silhouetteValues(const ICoreMatrix& points, const ICoreMatrix& assignment,
const std::string& metric, const double& exponent,
ICoreMatrix& out, std::string* whyNot = nullptr);
// ---- k-means, toolbox row T3.18 ---------------------------------------
//
// Only the arm that is DETERMINED: `kmeans(X, k, "Start", C0)` with the
// starting centres given as a matrix. Without them MATLAB picks its start
// from its own random stream (`"plus"`, k-means++) and no seed on this
// console reproduces the draw, so that form is refused rather than
// answered with a number nothing can compare (board rule Z11).
//
// Its five metric words are not five distance functions with a shared
// centre; each word carries its OWN centre, and that is the half of
// `kmeans` a reader is most likely to assume away:
//
// sqeuclidean squared distance, centre = the MEAN
// cityblock absolute distance, centre = the MEDIAN (the even case
// averages the two middle values, so a centre need not be
// one of the points)
// cosine 1 - the cosine of the angle, on rows normalized to unit
// length first; the centre is the UNNORMALIZED mean of the
// normalized rows, re-normalized only when it is measured
// against
// correlation the same, on rows centred on their own mean first
// hamming the fraction of coordinates that differ, on 0/1 data;
// centre = the MAJORITY, `2*sum > count`, so a tie goes to
// zero
enum class Centre { SquaredEuclidean, CityBlock, Cosine, Correlation, Hamming };
static bool centreMetricOf(const std::string& word, Centre& metric);
// `[idx, C, sumd, D] = kmeans(X, k, "Start", C0[, "Distance", d]
// [, "MaxIter", n])`.
//
// Two rules decide the answer where a reimplementation would guess:
// a point moves to another cluster only when the new distance is STRICTLY
// smaller (so a tie leaves it where it is), and `min` takes the FIRST of
// several equal distances (so a first assignment ties to the lower cluster
// number). An empty cluster is refilled with the single loneliest point,
// taken out of the cluster it was in -- MATLAB's `"singleton"` action,
// which is its default and the only one carried here.
static bool kMeans(const ICoreMatrix& points, const ICoreMatrix& start, const Centre& metric,
const size_t& maxIterations, ICoreMatrix& assignment, ICoreMatrix& centres,
ICoreMatrix& withinSums, ICoreMatrix& distances,
std::string* whyNot = nullptr);
// `[idx, C, sumd, D, midx] = kmedoids(X, k, "Start", C0[, "Distance", d])`
// -- the same question with the centre required to BE one of the points.
//
// ⚠ Its starting rows must be rows OF X, which `kmeans` does not require
// and MATLAB enforces with its own sentence: a medoid that is not a data
// point is not a medoid. `midx` is which rows they are.
//
// The swap phase of PAM, from the given start: every (medoid, non-medoid)
// exchange is costed, the best improving one is taken, and it repeats
// until none improves. That is deterministic given the start, which is why
// this arm is comparable and the unstarted one is not. `start` here is the
// metric's own distance -- squared Euclidean by default, which is what
// makes `sumd` a sum of SQUARED distances.
static bool kMedoids(const ICoreMatrix& points, const ICoreMatrix& start, const Centre& metric,
const size_t& maxIterations, ICoreMatrix& assignment,
ICoreMatrix& centres, ICoreMatrix& withinSums, ICoreMatrix& distances,
ICoreMatrix& medoidRows, std::string* whyNot = nullptr);
};
};ICoreCorrelationAnalysis.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreCorrelationAnalysis.h
ICoreCorrelationAnalysis#
ICoreCorrelationAnalysis.h:49 · class · 3 declaration(s)
The correlation block of toolbox row T3.16 -- corr, partialcorr, corrcov and tiedrank.
class ICoreCorrelationAnalysis {
public:
enum class Type { Pearson, Spearman, Kendall };
enum class Rows { All, Complete, Pairwise };
enum class Tail { Both, Right, Left };
static bool typeOf(const std::string& word, Type& type);
static bool rowsOf(const std::string& word, Rows& rows);
static bool tailOf(const std::string& word, Tail& tail);
// `[rho, pval] = corr(X[, Y])`. `y` may be empty for the one-argument
// form, which answers the symmetric matrix with exact ones on its
// diagonal (MATLAB puts them there rather than letting the arithmetic
// decide). `wantP` asks for the p-value and is refused for the two rank
// types, for the reason in the class comment.
static bool correlation(const ICoreMatrix& x, const ICoreMatrix* y, const Type& type,
const Rows& rows, const Tail& tail, const bool& wantP,
ICoreMatrix& rho, ICoreMatrix& pValue,
std::string* whyNot = nullptr);
// `[rho, pval] = partialcorr(X[, Y], Z)` -- the correlation of what is
// left of X (and Y) after regressing each column on `[1 Z]`. The p-value
// is the same t transform with `n - rank(Z) - 2` degrees of freedom.
static bool partialCorrelation(const ICoreMatrix& x, const ICoreMatrix* y,
const ICoreMatrix& z, const Type& type, const Tail& tail,
const bool& wantP, ICoreMatrix& rho, ICoreMatrix& pValue,
std::string* whyNot = nullptr);
// `[R, sigma] = corrcov(C)` -- the correlation matrix of a covariance
// matrix, and the standard deviations it was scaled by.
static bool correlationFromCovariance(const ICoreMatrix& c, ICoreMatrix& r,
ICoreMatrix& sigma, std::string* whyNot = nullptr);
// `[r, tieadj] = tiedrank(x)` -- ranks down each column, ties averaged.
// `tieadj` is the adjustment `sum((t^3 - t) / 2)` over the tie groups,
// which is what `signrank` and `ranksum` need and is 0 when nothing ties.
static bool tiedRank(const ICoreMatrix& x, ICoreMatrix& ranks, double& tieAdjustment,
std::string* whyNot = nullptr);
};
};ICoreDescriptiveStatistics.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreDescriptiveStatistics.h
ICoreDescriptiveStatistics#
ICoreDescriptiveStatistics.h:52 · class · 0 declaration(s)
The nine descriptive statistics of toolbox row T3.1 -- prctile, quantile, iqr, range, mad, zscore, skewness, kurtosis and moment.
class ICoreDescriptiveStatistics {
public:
// `prctile(x, p)` with p in PERCENT, and `quantile(x, p)` with p in [0, 1]
// -- one function because quantile.m is one line calling prctile with
// `100.*p`. `p` may be a vector, and the answer follows MATLAB's shape:
// a vector `x` answers p's shape, a matrix answers one row per percentile.
//
// `asQuantile` also carries quantile's own integer form: a scalar whole
// `p > 1` means "that many evenly spaced quantiles", `(1:p)/(p+1)`.
static ICoreMatrix percentile(const ICoreMatrix& x, const ICoreMatrix& p,
const int& dim, const bool& asQuantile,
std::string* whyNot = nullptr);
// `iqr(x)` -- the 75th percentile less the 25th, through the rule above,
// so it inherits the clamping and the (i-0.5)/n placement.
static ICoreMatrix interquartileRange(const ICoreMatrix& x, const int& dim,
std::string* whyNot = nullptr);
// `range(x)` -- max less min. Statistics Toolbox, though it is two core
// calls; the console keeps `max(A) - min(A)` as the core spelling.
static ICoreMatrix range(const ICoreMatrix& x, const int& dim,
std::string* whyNot = nullptr);
// `mad(x, flag)`. `useMedian` is MATLAB's `flag == 1`.
static ICoreMatrix absoluteDeviation(const ICoreMatrix& x, const bool& useMedian,
const int& dim, std::string* whyNot = nullptr);
// `[z, mu, sigma] = zscore(x, flag, dim)`. `population` is flag == 1 (N);
// the default is flag == 0 (N-1).
static bool standardScore(const ICoreMatrix& x, const bool& population, const int& dim,
ICoreMatrix& z, ICoreMatrix& mu, ICoreMatrix& sigma,
std::string* whyNot = nullptr);
// `skewness(x, flag)` and `kurtosis(x, flag)`. `biased` is flag == 1, and
// it is the DEFAULT. The unbiased correction is undefined below n = 3
// (skewness) and n = 4 (kurtosis) and answers NaN there, as MATLAB does.
static ICoreMatrix skewness(const ICoreMatrix& x, const bool& biased, const int& dim,
std::string* whyNot = nullptr);
static ICoreMatrix kurtosis(const ICoreMatrix& x, const bool& biased, const int& dim,
std::string* whyNot = nullptr);
// `moment(x, order)` -- the CENTRAL moment, always over the biased mean.
// Order 1 is exactly 0 by definition and MATLAB returns it without
// computing anything, which is why it is 0 and not 1e-17.
static ICoreMatrix centralMoment(const ICoreMatrix& x, const double& order, const int& dim,
std::string* whyNot = nullptr);
};
};ICoreDistanceMetrics.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreDistanceMetrics.h
ICoreDistanceMetrics#
ICoreDistanceMetrics.h:45 · class · 2 declaration(s)
The distance block of toolbox row T3.19 -- pdist, pdist2, squareform, mahal and knnsearch, over twelve of MATLAB's distance metrics.
class ICoreDistanceMetrics {
public:
// MATLAB's twelve metric words. `Jaccard` and `Hamming` read their inputs
// as patterns rather than numbers, which is why they are here rather than
// refused: MATLAB defines both on any numeric data.
enum class Metric {
Euclidean, SquaredEuclidean, StandardizedEuclidean, CityBlock, Minkowski,
Chebychev, Cosine, Correlation, Hamming, Jaccard, Spearman, Mahalanobis
};
// The metric word, MATLAB's spellings (including "chebychev", which is
// MATLAB's spelling and not the more common "chebyshev").
static bool metricOf(const std::string& word, Metric& metric);
// `pdist(X[, metric[, arg]])` -- the CONDENSED 1 x n(n-1)/2 row.
// `exponent` is Minkowski's p and is ignored by every other metric.
static ICoreMatrix pairwise(const ICoreMatrix& x, const Metric& metric,
const double& exponent, std::string* whyNot = nullptr);
// `pdist2(X, Y[, metric[, arg]])` -- the FULL nx x ny matrix. The scale
// (standardized Euclidean) and the covariance (Mahalanobis) come from X
// alone, which is MATLAB's default and is measurable rather than
// documented: pdist2(X, Y, "seuclidean") is pdist's own weights.
static ICoreMatrix pairwiseBetween(const ICoreMatrix& x, const ICoreMatrix& y,
const Metric& metric, const double& exponent,
std::string* whyNot = nullptr);
// `squareform(v)` and `squareform(D)`, the direction taken from the shape
// exactly as squareform.m takes it: a VECTOR becomes the symmetric matrix,
// a matrix becomes the vector, and a matrix with a non-zero diagonal is an
// error rather than a matrix with its diagonal ignored.
static ICoreMatrix squareForm(const ICoreMatrix& in, std::string* whyNot = nullptr);
// `mahal(Y, X)` -- the SQUARED Mahalanobis distance of each row of Y from
// the sample X, one per row of Y. Not the same number as
// `pdist2(Y, mean(X), "mahalanobis")`, which is its square root.
static ICoreMatrix mahalanobis(const ICoreMatrix& y, const ICoreMatrix& x,
std::string* whyNot = nullptr);
// `[idx, D] = knnsearch(X, Y[, "K", k][, "Distance", metric])` -- the k
// nearest rows of X to each row of Y, nearest first, as a ny x k index
// matrix (1-based, MATLAB's) and the matching distances.
//
// TIES GO TO THE SMALLER INDEX. Measured against R2026a rather than
// assumed: knnsearch([1 0; 0 1; 5 5], [0 0], "K", 2) answers [1 2], and
// MATLAB's own documentation only promises "one of them".
static bool nearestNeighbours(const ICoreMatrix& x, const ICoreMatrix& y, const size_t& k,
const Metric& metric, const double& exponent,
ICoreMatrix& indices, ICoreMatrix& distances,
std::string* whyNot = nullptr);
};
};ICoreDistributions.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreDistributions.h
ICoreDistributions#
ICoreDistributions.h:51 · class · nested FamilyShape · 1 declaration(s)
The normal family and the t / chi-square / F families (toolbox rows T3.3 and T3.5): the distribution FUNCTIONS, which take numbers and answer numbers.
class ICoreDistributions {
public:
// Whether an answer is the lower tail or MATLAB's "upper" one.
enum class Tail { Lower, Upper };
// ---- the normal family (T3.3) ----------------------------------------
//
// `x`, `mu` and `sigma` broadcast the way MATLAB's `distchck` broadcasts:
// each is either a scalar or the common shape. A `sigma` at or below zero
// answers NaN rather than refusing, which is MATLAB's choice throughout
// this family.
static bool normalPdf(const ICoreMatrix& x, const ICoreMatrix& mean,
const ICoreMatrix& deviation, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool normalCdf(const ICoreMatrix& x, const ICoreMatrix& mean,
const ICoreMatrix& deviation, const Tail& tail, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool normalInv(const ICoreMatrix& p, const ICoreMatrix& mean,
const ICoreMatrix& deviation, ICoreMatrix& out,
std::string* whyNot = nullptr);
// MATLAB's [m, v] = normstat(mu, sigma): the mean and the VARIANCE, so v
// is sigma^2 and not sigma.
static bool normalStat(const ICoreMatrix& mean, const ICoreMatrix& deviation,
ICoreMatrix& outMean, ICoreMatrix& outVariance,
std::string* whyNot = nullptr);
// MATLAB's nlogL = normlike([mu sigma], data): the NEGATIVE log
// likelihood, so a better fit is a smaller number.
static bool normalLike(const ICoreMatrix& parameters, const ICoreMatrix& data,
double& negativeLogLikelihood, std::string* whyNot = nullptr);
// MATLAB's [muhat, sigmahat, muci, sigmaci] = normfit(x[, alpha]).
//
// `sigmahat` is the UNBIASED deviation (divided by n - 1), and the two
// intervals are the exact ones rather than a normal approximation: the
// mean's is a t interval and the deviation's a chi-square interval, which
// is why this row and T3.5 were claimed together -- `normfit` cannot be
// finished without `tinv` and `chi2inv`.
static bool normalFit(const ICoreMatrix& x, const double& alpha, double& meanEstimate,
double& deviationEstimate, ICoreMatrix& meanInterval,
ICoreMatrix& deviationInterval, std::string* whyNot = nullptr);
// ---- the t, chi-square and F families (T3.5) -------------------------
static bool studentPdf(const ICoreMatrix& x, const ICoreMatrix& degrees, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool studentCdf(const ICoreMatrix& x, const ICoreMatrix& degrees, const Tail& tail,
ICoreMatrix& out, std::string* whyNot = nullptr);
static bool studentInv(const ICoreMatrix& p, const ICoreMatrix& degrees, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool chiSquarePdf(const ICoreMatrix& x, const ICoreMatrix& degrees, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool chiSquareCdf(const ICoreMatrix& x, const ICoreMatrix& degrees, const Tail& tail,
ICoreMatrix& out, std::string* whyNot = nullptr);
static bool chiSquareInv(const ICoreMatrix& p, const ICoreMatrix& degrees, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool fisherPdf(const ICoreMatrix& x, const ICoreMatrix& first,
const ICoreMatrix& second, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool fisherCdf(const ICoreMatrix& x, const ICoreMatrix& first,
const ICoreMatrix& second, const Tail& tail, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool fisherInv(const ICoreMatrix& p, const ICoreMatrix& first,
const ICoreMatrix& second, ICoreMatrix& out,
std::string* whyNot = nullptr);
// ---- the four samplers (T3.3's normrnd, T3.5's trnd/chi2rnd/frnd) ----
//
// ⚠ **These four cannot be parity-compared at all, and it is not a
// tolerance question.** MATLAB's `randn` is a Mersenne Twister stream and
// this console's is not, so the same seed does not produce the same draw
// and no band makes it. They cross with `Compare::Skip` (C0.9) and are
// pinned by regress cases on a DERIVED property instead -- a count, a
// shape, a sign -- which is the shape `rand`'s own cases already take
// (session.itest: `numel(rand(2, 3))`, `nnz(floor(rand(5, 5)))`).
//
// `rows` and `columns` of 0 ask for MATLAB's parameter-shaped answer: a
// draw per entry of the broadcast parameters, which for scalars is one.
static bool normalRandom(const ICoreMatrix& mean, const ICoreMatrix& deviation,
const size_t& rows, const size_t& columns, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool studentRandom(const ICoreMatrix& degrees, const size_t& rows,
const size_t& columns, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool chiSquareRandom(const ICoreMatrix& degrees, const size_t& rows,
const size_t& columns, ICoreMatrix& out,
std::string* whyNot = nullptr);
static bool fisherRandom(const ICoreMatrix& first, const ICoreMatrix& second,
const size_t& rows, const size_t& columns, ICoreMatrix& out,
std::string* whyNot = nullptr);
// ---- the two inverses base MATLAB has and this console did not --------
//
// `betaincinv(p, a, b)` and `gammaincinv(p, a)`, solved by Newton on
// `betainc`/`gammainc`. They are PUBLIC because three of the six inverses
// above are one call to them plus a change of variable, and because the
// gamma and beta distribution families (toolbox row T3.6) will want the
// same two functions -- one implementation, so a tolerance cannot drift
// between the rows that share it.
static double incompleteBetaInverse(const double& p, const double& a, const double& b,
const Tail& tail);
static double incompleteGammaInverse(const double& p, const double& a,
const Tail& tail = Tail::Lower);
// ==== mb143 T3.4 / T3.6 begin ========================================
//
// ---- the remaining univariate families (toolbox rows T3.4 and T3.6) --
//
// Fifteen families, and they are ONE entry point each rather than ninety
// methods, because ninety methods would be ninety copies of the same four
// lines. What differs between `unifpdf` and `hygepdf` is a scalar formula
// and how many parameters come after x; what does NOT differ is MATLAB's
// `distchck` broadcast, the NaN-rather-than-error convention for a
// parameter out of range, the trailing `"upper"` on every cdf, and the
// `[m, n]` size tail on every sampler. Those four live here once, so a
// family cannot quietly get a different one -- the same argument that made
// `incompleteBetaInverse` public for this row rather than copied into it.
//
// The families the T3.5 header already named as this row's callers are
// exactly the ones that appear below as one line: `gamcdf(x, a, b)` IS
// `gammainc(x/b, a)` and `gaminv` IS `gammaincinv(p, a) * b`, `betacdf`
// and `betainv` are `betainc`/`betaincinv` unchanged.
enum class Family {
// T3.4
Uniform, Exponential, Poisson, Binomial,
// T3.6
Gamma, Beta, Lognormal, Weibull, Rayleigh, ExtremeValue,
GeneralizedExtremeValue, GeneralizedPareto, NegativeBinomial, Geometric,
Hypergeometric,
};
// How many parameters follow x, how many of them MATLAB insists on for a
// `*pdf`/`*cdf`/`*inv`/`*rnd`, and how many it insists on for the `*stat`
// pair -- which is NOT always the same number: `gpstat(k, sigma)` takes
// two of three and `gevstat(k, sigma, mu)` takes all three.
struct FamilyShape {
size_t parameters = 0;
size_t required = 0;
size_t requiredForMoments = 0;
// A count-valued family: its quantile is a whole number found by
// walking the cdf, not a closed form.
bool discrete = false;
};
static FamilyShape shapeOf(const Family& family);
// `arguments` is x (or p) followed by the family's parameters; a parameter
// left off takes MATLAB's default where it has one. Every argument is a
// scalar or the common shape.
static bool density(const Family& family, const std::vector<ICoreMatrix>& arguments,
ICoreMatrix& out, std::string* whyNot = nullptr);
static bool cumulative(const Family& family, const std::vector<ICoreMatrix>& arguments,
const Tail& tail, ICoreMatrix& out, std::string* whyNot = nullptr);
static bool quantile(const Family& family, const std::vector<ICoreMatrix>& arguments,
ICoreMatrix& out, std::string* whyNot = nullptr);
// MATLAB's [m, v] = <family>stat(...): the mean and the VARIANCE.
static bool moments(const Family& family, const std::vector<ICoreMatrix>& parameters,
ICoreMatrix& outMean, ICoreMatrix& outVariance,
std::string* whyNot = nullptr);
// ⚠ Every sampler here crosses the parity suite with `Compare::Skip` for
// `normrnd`'s reason -- MATLAB's `rand`/`randn`/`randg` are one Mersenne
// Twister stream and this console's are not -- so what these reproduce is
// the DISTRIBUTION and the shape, never the draw (Z11).
static bool sample(const Family& family, const std::vector<ICoreMatrix>& parameters,
const size_t& rows, const size_t& columns, ICoreMatrix& out,
std::string* whyNot = nullptr);
// MATLAB's `<family>fit(x[, alpha])`, and `binofit(x, n[, alpha])` -- the
// only one with a second DATA argument, which is why `trials` is a
// parameter here rather than folded into `data`.
//
// `estimate` is MATLAB's `parmhat` and `interval` its `parmci`, each in
// MATLAB's own shape: `expfit` answers a number and a 2 x 1, `lognfit` a
// 1 x 2 and a 2 x 2, `binofit` a number and a 1 x 2, and `unifit` a 1 x 2
// whose interval columns are its third and fourth outputs.
//
// The four families whose fit is `fminsearch` on the log-likelihood
// (`beta`, `nbin`, `gev`, `gp`) are NOT here and are refused by name: a
// simplex that stops at TolX 1e-6 reports the point it stopped at, and
// that is determined by the PATH rather than by an equation (C8.22). The
// four that are here are a closed form or a bracketed root, which is not.
static bool fit(const Family& family, const ICoreMatrix& data, const ICoreMatrix& trials,
const double& alpha, ICoreMatrix& estimate, ICoreMatrix& interval,
std::string* whyNot = nullptr);
// ==== mb143 T3.4 / T3.6 end ==========================================
// ==== mb143 T3.9 begin ===============================================
//
// ---- `mle` for a NAMED distribution (toolbox row T3.9) ---------------
//
// ⚠ **`mle` is not an estimator.** `mle.m` is a switch over eighteen
// distribution names that calls the `*fit` functions above -- `binofit`,
// `unifit`, `expfit`, `poissfit`, `normfit`, `lognfit`, `raylfit`,
// `gamfit`, `evfit`, `wblfit` -- plus four closed forms of its own, and
// routes everything else through `fitdist` to a distribution OBJECT. So
// this row is a dispatch, and the only arithmetic it adds is those four
// and one correction.
//
// ⚠ **The correction is the trap, and it is a real one**:
// `mle(x, "Distribution", "normal")` is NOT `normfit(x)`. `fitdist`
// answers the UNBIASED deviation (divided by n - 1) and `mle` multiplies
// it by `sqrt((n - 1)/n)` to get the maximum-likelihood one, so the two
// functions disagree on the second parameter by 6.5% at n = 8 and agree
// on the first exactly. Measured on R2026a: `normfit` 0.86147713674992976,
// `mle` 0.80583807306430988. **The INTERVAL is not corrected** -- `mle`
// hands back `normfit`'s, built from the unbiased estimate -- so a reader
// who "fixes" the interval to match the estimate breaks it. Lognormal
// does the same thing on the same line of `mle.m`.
//
// Three of the names here have no `*pdf` on this console at all, because
// no row asked for one and MATLAB has no `bernfit`, `unidfit` or `hnfit`
// either: `mle` is the only way to fit them, which is why they are an
// enumerator of their own rather than three more `Family` members whose
// density nothing would answer.
enum class MleDistribution {
Normal, Lognormal, Exponential, Poisson, Rayleigh, Gamma, Weibull,
ExtremeValue, Uniform, Geometric, DiscreteUniform, Bernoulli, Binomial,
HalfNormal,
};
// `trials` is `mle`'s `"ntrials"` (binomial only) and `location` its
// `"mu"` (half normal only); both are ignored elsewhere, as MATLAB
// ignores them with a warning. `estimate` and `interval` are MATLAB's
// `phat` and `pci` in MATLAB's own shapes -- and note that a
// one-parameter family's `pci` is a 2 x 1 here where `binofit`'s is a
// 1 x 2, because `mle.m` transposes it.
static bool maximumLikelihood(const MleDistribution& distribution, const ICoreMatrix& data,
const double& alpha, const double& trials,
const double& location, ICoreMatrix& estimate,
ICoreMatrix& interval, std::string* whyNot = nullptr);
// ==== mb143 T3.9 end =================================================
};
};ICoreEmpiricalDistributions.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreEmpiricalDistributions.h
ICoreEmpiricalDistributions#
ICoreEmpiricalDistributions.h:47 · class · nested KernelRequest · 4 declaration(s)
The two empirical distribution estimates of toolbox row T3.24 -- MATLAB's ecdf and ksdensity.
class ICoreEmpiricalDistributions {
public:
// MATLAB's three `ecdf` functions. `"chf"` and `"cumulative hazard"` are
// spellings of the third.
enum class Empirical { Cdf, Survivor, CumulativeHazard };
// The function word, MATLAB's own list including its abbreviations.
static bool empiricalOf(const std::string& word, Empirical& which);
// `[f, x, flo, fup] = ecdf(y, ...)`. `censoring` and `frequency` may be
// empty matrices, which is "none given". `alpha` is the 0.05 of MATLAB's
// default and only matters when `wantBounds` is true; the bounds are
// Greenwood's, with a NaN in the first row because the starting point is
// not an estimate.
//
// Left-censored, interval-censored and doubly censored data are refused:
// each is a different estimator in MATLAB (Ware-Demets, Turnbull, and an
// EM iteration), not an option of this one.
static bool empiricalCdf(const ICoreMatrix& y, const ICoreMatrix& censoring,
const ICoreMatrix& frequency, const Empirical& which,
const double& alpha, const bool& wantBounds,
ICoreMatrix& f, ICoreMatrix& x,
ICoreMatrix& lower, ICoreMatrix& upper,
std::string* whyNot = nullptr);
// MATLAB's four built-in smoothing kernels. Each is scaled to unit
// variance, which is what makes the bandwidth mean the same thing across
// them -- and why each has its own cutoff.
enum class Kernel { Normal, Box, Triangle, Epanechnikov };
// What the estimate answers. `Icdf` is deliberately absent (see above).
enum class Density { Pdf, Cdf, Survivor, CumulativeHazard };
static bool kernelOf(const std::string& word, Kernel& kernel);
static bool densityOf(const std::string& word, Density& which);
// Everything a `ksdensity` call carries besides its data. `bandwidth` at
// or below zero means "estimate it", which is MATLAB's empty default;
// `lower`/`upper` are the support, `-Inf`/`Inf` for the unbounded default.
struct KernelRequest {
Kernel kernel = Kernel::Normal;
Density function = Density::Pdf;
double bandwidth = 0.0;
double lower = 0.0;
double upper = 0.0;
size_t points = 100;
bool pointsGiven = false;
ICoreMatrix at;
};
// A request with MATLAB's defaults already in it: normal kernel, pdf,
// estimated bandwidth, unbounded support, 100 points.
static KernelRequest defaultRequest();
// `[f, xi, u] = ksdensity(y, ...)`. `at` comes back holding the points the
// estimate was evaluated at -- the caller's own when it gave some, and the
// default grid otherwise -- and `bandwidth` the one that was used, which
// is the third output MATLAB answers.
static bool kernelDensity(const ICoreMatrix& y, const KernelRequest& request,
ICoreMatrix& f, ICoreMatrix& at, double& bandwidth,
std::string* whyNot = nullptr);
};
};ICoreGroupComparisons.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreGroupComparisons.h
ICoreGroupComparisons#
ICoreGroupComparisons.h:39 · class · 0 declaration(s)
The nonparametric tests and the analysis of variance of toolbox row T3.12 -- signtest, signrank, ranksum, anova1, anova2, kruskalwallis and friedman.
class ICoreGroupComparisons {
public:
using Tail = ICoreHypothesisTests::Tail;
// `[p, h] = signtest(x[, y][, "Alpha", a][, "Tail", t])` -- are the
// differences centred on zero? Only their SIGNS are used, so it is the
// weakest of the three and the one that survives any monotone rescaling.
//
// Pairs that differ by exactly zero are DROPPED, and the binomial the
// exact arm uses is over what is left.
static bool signTest(const ICoreMatrix& x, const ICoreMatrix& y, const bool& paired,
const double& alpha, const Tail& tail, double& probability,
double& decision, std::string* whyNot = nullptr);
// `[p, h] = signrank(x[, y][, "Alpha", a][, "Tail", t])` -- the same
// question using the RANKS of the differences as well as their signs.
//
// ⚠ Its zero test is not `== 0`. `signrank.m` drops a pair whose
// difference is within `eps(x) + eps(y)` of zero, which is a wider net
// than equality and is what stops two numbers that differ in their last
// bit from contributing a rank.
static bool signedRankTest(const ICoreMatrix& x, const ICoreMatrix& y, const bool& paired,
const double& alpha, const Tail& tail, double& probability,
double& decision, std::string* whyNot = nullptr);
// `[p, h] = ranksum(x, y[, "Alpha", a][, "Tail", t])` -- the two-sample
// rank sum (Mann-Whitney U), on samples that need not be the same length.
//
// The statistic is the rank sum of the SMALLER sample, so which of the two
// arguments that is decides what a one-sided test means; MATLAB tracks it
// and flips the tail to match, and so does this.
static bool rankSumTest(const ICoreMatrix& x, const ICoreMatrix& y, const double& alpha,
const Tail& tail, double& probability, double& decision,
std::string* whyNot = nullptr);
// `p = anova1(y[, group])` and `p = kruskalwallis(y[, group])` -- one-way,
// over the COLUMNS of a matrix or over a vector split by a grouping
// vector. `kruskalwallis` is `anova1` on the tied ranks with a
// chi-square in place of the F, which is how `kruskalwallis.m` is written:
// one line, forwarding to `anova1`.
static bool oneWay(const ICoreMatrix& x, const ICoreMatrix& group, const bool& grouped,
const bool& byRanks, double& probability, double& statistic,
double& firstDegrees, double& secondDegrees,
std::string* whyNot = nullptr);
// `p = anova2(y[, reps])` -- two-way, the rows being one factor and the
// columns the other. Answers MATLAB's row of p-values: the column effect,
// the row effect, and with more than one replicate the interaction too.
static bool twoWay(const ICoreMatrix& x, const size_t& replicates, ICoreMatrix& probabilities,
std::string* whyNot = nullptr);
// `p = friedman(y[, reps])` -- the rank analogue of `anova2`, ranking
// WITHIN each block before the analysis and reading a chi-square off the
// column sum of squares.
static bool friedmanTest(const ICoreMatrix& x, const size_t& replicates, double& probability,
double& statistic, std::string* whyNot = nullptr);
};
};ICoreHypothesisTests.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreHypothesisTests.h
ICoreHypothesisTests#
ICoreHypothesisTests.h:42 · class · 2 declaration(s)
The variance and goodness-of-fit tests of toolbox row T3.11 -- vartest, vartest2, chi2gof, kstest, kstest2, lillietest, jbtest, adtest and runstest.
class ICoreHypothesisTests {
public:
// ⚠ THE TAIL WORDS ARE NOT THE SAME ACROSS THIS ROW, and neither is the
// CLASS of the decision. `vartest`, `vartest2` and `runstest` take
// `"both"`, `"left"` and `"right"` and answer a DOUBLE `h`; `kstest` and
// `kstest2` take `"unequal"`, `"smaller"` and `"larger"` for the same
// three tails and answer a LOGICAL one, as does `adtest`. Both
// inconsistencies are MATLAB's and both are carried: a console that took
// `"both"` for `kstest` would accept a script MATLAB refuses, and one
// that answered a double there would report a different `class()` for the
// same call. `tailOf` below is the first vocabulary; the second lives at
// the call site, because it belongs to two names rather than to the row.
enum class Tail { Left, Both, Right };
static bool tailOf(const std::string& word, Tail& tail);
// The three distribution families `lillietest` and `adtest` share. MATLAB
// reaches the lognormal and the Weibull by taking logs and testing the
// normal and the extreme value, which is what this does.
enum class Family { Normal, Exponential, ExtremeValue, Lognormal, Weibull };
static bool familyOf(const std::string& word, Family& family);
// `[h, p, ci] = vartest(x, v[, "Alpha", a][, "Tail", t])` -- is the
// variance of x equal to v? The statistic is `sum((x - mean(x))^2) / v`
// against a chi-square with n-1 degrees of freedom.
static bool varianceTest(const ICoreMatrix& x, const double& variance, const double& alpha,
const Tail& tail, double& decision, double& probability,
ICoreMatrix& interval, std::string* whyNot = nullptr);
// `[h, p, ci] = vartest2(x, y[, "Alpha", a][, "Tail", t])` -- do two
// samples have the same variance? The ratio of the two variances against
// an F.
static bool varianceTest2(const ICoreMatrix& x, const ICoreMatrix& y, const double& alpha,
const Tail& tail, double& decision, double& probability,
ICoreMatrix& interval, std::string* whyNot = nullptr);
// `[h, p, stat, df] = chi2gof(x[, "NBins", n][, "Emin", m][, "NParams", k]
// [, "Alpha", a])`.
//
// ⚠ THE POOLING IS THE TEST. MATLAB bins into ten equal-width bins by
// default, then POOLS neighbouring bins from whichever end has the
// smaller expected count until every remaining bin expects at least
// `emin` (5 by default) -- so the degrees of freedom depend on the data,
// not only on the bin count. A `chi2gof` that skips the pooling answers a
// different number with the same formula.
//
// The default null is a normal fitted to the data, so two parameters are
// estimated and `df = bins - 1 - 2`.
static bool chiSquareGoodnessOfFit(const ICoreMatrix& x, const size_t& bins,
const double& minimumExpected, const double& parameters,
const double& alpha, double& decision, double& probability,
double& statistic, double& degrees,
std::string* whyNot = nullptr);
// `[h, p, ksstat, cv] = kstest(x[, "Alpha", a][, "Tail", t])` against the
// STANDARD normal -- which is MATLAB's default null and catches people
// out: `kstest` does not standardize x for you, and `kstest(x)` on data
// with any other mean or spread rejects for that reason alone.
//
// The p-value is exact rather than asymptotic below a threshold MATLAB
// sets on `n * D^2`: Marsaglia, Tsang and Wang's method, which raises a
// (2k-1)-square Toeplitz matrix to the n-th power. Above the threshold it
// is a three-term asymptotic form. Both arms are here, because the
// threshold is part of the answer.
//
// ⚠ The CRITICAL VALUE has three arms of its own, chosen on the ONE-SIDED
// level: outside [0.005, 0.10] MATLAB answers NaN rather than
// extrapolating; at or under twenty observations it reads an exact table
// with a NOT-A-KNOT SPLINE -- a different interpolant from the `pchip`
// the three tabulated tests use, on the same kind of table -- and above
// twenty it is a three-term fit capped at `1 - alpha`.
static bool kolmogorovSmirnov(const ICoreMatrix& x, const double& alpha, const Tail& tail,
double& decision, double& probability, double& statistic,
double& critical, std::string* whyNot = nullptr);
// `[h, p, ks2stat] = kstest2(x, y[, "Alpha", a][, "Tail", t])` -- the
// two-sample test, whose p-value is the asymptotic Kolmogorov series with
// Stephens' small-sample correction on lambda.
static bool kolmogorovSmirnov2(const ICoreMatrix& x, const ICoreMatrix& y, const double& alpha,
const Tail& tail, double& decision, double& probability,
double& statistic, std::string* whyNot = nullptr);
// `[h, p, kstat, critval] = lillietest(x[, "Distribution", d][, "Alpha", a])`
// -- the Kolmogorov-Smirnov statistic against a distribution whose
// parameters were ESTIMATED FROM THE SAME DATA, which is why its null
// distribution is tabulated rather than computed.
static bool lilliefors(const ICoreMatrix& x, const Family& family, const double& alpha,
double& decision, double& probability, double& statistic,
double& critical, std::string* whyNot = nullptr);
// `[h, p, adstat, cv] = adtest(x[, "Distribution", d][, "Alpha", a])` --
// the Anderson-Darling statistic, which weights the tails where the
// Kolmogorov-Smirnov statistic weights the middle. Same reason for the
// table, and its table is interpolated in LOG alpha rather than in alpha.
static bool andersonDarling(const ICoreMatrix& x, const Family& family, const double& alpha,
double& decision, double& probability, double& statistic,
double& critical, std::string* whyNot = nullptr);
// `[h, p, jbstat, critval] = jbtest(x[, alpha])` -- skewness and kurtosis
// together, against a table indexed by sample size AND by alpha, which
// MATLAB interpolates in `1/n` first and then in alpha.
static bool jarqueBera(const ICoreMatrix& x, const double& alpha, double& decision,
double& probability, double& statistic, double& critical,
std::string* whyNot = nullptr);
// `[h, p] = runstest(x[, v][, "Alpha", a][, "Tail", t])` -- are the values
// above and below `v` in random order?
//
// ⚠ THE DEFAULT `v` IS THE MEAN, not the median. `runstest.m` says
// `v = mean(x)` and a reader who assumes the median gets a different
// sequence of ones and zeros and a different answer. Values exactly equal
// to `v` are DROPPED before anything is counted.
//
// The p-value is exact -- a closed form in binomial coefficients over the
// number of runs -- and not the normal approximation, which MATLAB keeps
// only for the up-and-down variant this console refuses.
static bool runsTest(const ICoreMatrix& x, const double& threshold, const double& alpha,
const Tail& tail, double& decision, double& probability,
double& runs, double& normalScore, std::string* whyNot = nullptr);
// ---- the LOCATION tests (toolbox row T3.10) --------------------------
//
// `ttest`, `ttest2` and `ztest`: three tests of a mean, and every number
// they answer is a composition over cdfs T3.5 and T3.3 already landed --
// nothing here is estimated or tabulated. What the row is really about is
// three conventions that no signature states:
//
// ⚠ `ttest(x, y)` IS PAIRED. Two samples handed to `ttest` are the
// SAME subjects measured twice, and the test runs on x - y with n - 1
// degrees of freedom; `ttest2(x, y)` is the unpaired test with
// n1 + n2 - 2. The two answer different p-values for the same pair of
// vectors and neither call says which it is.
//
// ⚠ A ONE-SIDED interval is HALF-INFINITE. MATLAB's right-tailed ci is
// [lower, Inf] and its left-tailed one [-Inf, upper] -- not the
// two-sided interval with one end reported.
//
// ⚠ 'Vartype', 'unequal' changes the DEGREES OF FREEDOM and not just
// the standard error: Welch-Satterthwaite gives a fractional df
// (14.21215207001414 on this row's own fixture), and a test that
// kept n1 + n2 - 2 there answers a p-value that is wrong in the
// third digit and plausible in every other way.
//
// `stats` -- MATLAB's fourth output, a struct of tstat/df/sd -- is a kind
// this console does not have (T0.6) and is refused by name at the call
// site rather than flattened into a vector whose fields nobody can tell
// apart.
//
// `paired` selects ttest's two-sample reading; `pooled` selects ttest2's
// (false is Welch). `deviation` above zero makes it `ztest` -- the
// normal test with a KNOWN standard deviation, where the interval is
// MATLAB's z-interval and the p-value the normal tail.
static bool locationTest(const ICoreMatrix& x, const ICoreMatrix& y, const double& mean,
const double& alpha, const Tail& tail, const bool& paired,
const bool& pooled, const double& deviation, double& decision,
double& probability, ICoreMatrix& interval, double& statistic,
double& degreesOfFreedom, double& standardDeviation,
std::string* whyNot = nullptr);
};
};ICoreLinearRegression.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreLinearRegression.h
ICoreLinearRegression#
ICoreLinearRegression.h:50 · class · 0 declaration(s)
regress -- toolbox row T3.13, and all five of its outputs.
class ICoreLinearRegression {
public:
// `[b, bint, r, rint, stats] = regress(y, X[, alpha])`. Rows where `y` or
// any column of `X` is missing are dropped before anything is computed and
// NaN is put back into `r` and `rint` afterwards, as regress.m does.
static bool regress(const ICoreMatrix& y, const ICoreMatrix& x, const double& alpha,
ICoreMatrix& coefficients, ICoreMatrix& coefficientInterval,
ICoreMatrix& residuals, ICoreMatrix& residualInterval,
ICoreMatrix& statistics, std::string* whyNot = nullptr);
};
};ICoreMultivariateAnalysis.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreMultivariateAnalysis.h
ICoreMultivariateAnalysis#
ICoreMultivariateAnalysis.h:34 · class · 1 declaration(s)
The multivariate block of the Statistics and Machine Learning Toolbox that is NOT principal components -- toolbox row T3.26: canoncorr, cmdscale, procrustes and rotatefactors.
class ICoreMultivariateAnalysis {
public:
// `[A, B, r, U, V] = canoncorr(X, Y)` -- the canonical correlations of two
// sets of columns measured on the same n rows. Both are centred first;
// `A` and `B` are the weights that carry them into the canonical
// variables `U = Xc*A` and `V = Yc*B`, and `r(k)` is the correlation
// between column k of each.
//
// A rank-deficient X or Y is NOT refused: `canoncorr.m` warns, drops the
// deficient columns from the factorization and puts explicit zero rows
// back into A and B at the pivot positions, and this does the same, so a
// constant column answers a zero weight rather than an error.
static bool canonicalCorrelation(const ICoreMatrix& x, const ICoreMatrix& y,
ICoreMatrix& a, ICoreMatrix& b, ICoreMatrix& r,
ICoreMatrix& u, ICoreMatrix& v,
std::string* whyNot = nullptr);
// `[Y, e] = cmdscale(D)` and `cmdscale(D, p)` -- classical (metric)
// multidimensional scaling: the point configuration whose interpoint
// distances reproduce D.
//
// `d` may be a full square dissimilarity matrix, a SIMILARITY matrix (unit
// diagonal and every entry below 1, which cmdscale.m turns into
// `sqrt(1 - D)`), or the `pdist` row of the lower triangle, which it
// squareforms. `dimensions` is 0 for "all that survive"; a positive p asks
// for exactly p columns.
//
// ⚠ `e` IS THE WHOLE EIGENVALUE SPECTRUM AND `Y` IS NOT ALL OF IT: the
// columns kept are those with `e > max(abs(e)) * eps^(3/4)`, so a set of n
// points in 3 dimensions answers an n x 3 Y beside an n x n `e` whose tail
// is rounding noise. Reading `size(Y, 2)` as `numel(e)` is the trap.
static bool classicalScaling(const ICoreMatrix& d, const size_t& dimensions,
ICoreMatrix& configuration, ICoreMatrix& eigenvalues,
std::string* whyNot = nullptr);
// `procrustes.m`'s 'Reflection': 'best' lets the data decide, and the two
// explicit values force one in or out.
enum class Reflection { Best, Forced, Forbidden };
// `[d, Z, transform] = procrustes(X, Y)` -- the rotation, scaling and
// translation of Y that lands it closest to X, and the standardized
// residual `d` that is left.
//
// `transform` is a struct in MATLAB and this console has no struct kind
// (T0.6), so the CONSOLE's `procrustes` stops at `[d, Z]` and refuses the
// third output by name -- inventing a three-output spelling MATLAB does
// not have would be the silent divergence C0.1 rejects. This function
// still answers the three fields, because the parity helper and the
// console's own tests need them: `rotation` is `transform.T`, `scale` is
// `transform.b` and `translation` is the one ROW of `transform.c` (MATLAB
// repeats that row n times).
//
// ⚠ Y MAY HAVE FEWER COLUMNS THAN X and is zero-padded up to X's
// dimension, which is how a 2-D shape is fitted into a 3-D one; `T` then
// comes back with Y's own number of ROWS, not X's.
static bool procrustesAnalysis(const ICoreMatrix& x, const ICoreMatrix& y,
const bool& scaling, const Reflection& reflection,
double& dissimilarity, ICoreMatrix& transformed,
ICoreMatrix& rotation, double& scale,
ICoreMatrix& translation, std::string* whyNot = nullptr);
// `rotatefactors.m`'s methods, less the two this row does not land.
// Orthomax is the family: varimax is gamma = 1, quartimax gamma = 0,
// equamax gamma = m/2 and parsimax gamma = d*(m-1)/(d+m-2), all four
// reaching the same iteration through a different coefficient.
enum class Rotation { Orthomax, Varimax, Quartimax, Equamax, Parsimax,
Procrustes, Promax };
// `[B, T] = rotatefactors(A, ...)`. `coefficient` is 'Coeff' (used only by
// Orthomax), `power` is 'Power' (Promax, default 4), `normalize` is
// 'Normalize' (default on -- Kaiser's row normalisation), `tolerance` is
// 'Reltol' (default sqrt(eps)) and `iterations` is 'Maxit' (default 250).
// `target` and `oblique` are the Procrustes method's 'Target' and 'Type';
// `target` is ignored by every other method.
//
// ⚠ ORTHOMAX CAN ASK FOR A RANDOM START AND THIS REFUSES INSTEAD. When the
// identity rotation already satisfies the convergence test at the first
// step -- a loading matrix that is ALREADY rotated -- `rotatefactors.m`
// restarts from `qr(randn(m, m))`, so its answer is a draw from MATLAB's
// own Mersenne Twister and reproducible on no other stream (Z11). That
// branch refuses by name; every other input takes the deterministic path
// from `T = eye(m)` and agrees with MATLAB entry for entry.
static bool rotateFactors(const ICoreMatrix& loadings, const Rotation& method,
const double& coefficient, const double& power,
const bool& normalize, const double& tolerance,
const size_t& iterations, const ICoreMatrix& target,
const bool& oblique, ICoreMatrix& rotated,
ICoreMatrix& rotation, std::string* whyNot = nullptr);
// The method word, as `rotatefactors(A, 'Method', word)` spells it.
// 'equimax' and 'equamax' are the same rotation, which rotatefactors.m
// accepts under both spellings and under the prefix 'equ'.
static bool rotationOf(const std::string& word, Rotation& out);
};
};ICoreMultivariateDensities.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreMultivariateDensities.h
ICoreMultivariateDensities#
ICoreMultivariateDensities.h:36 · class · 0 declaration(s)
The multivariate normal and t densities of toolbox row T3.8 -- mvnpdf, mvncdf, mvtpdf and the four samplers that go with them.
class ICoreMultivariateDensities {
public:
// `mvnpdf(X[, Mu[, Sigma]])` -- one density per ROW of X, as a column.
//
// `mu` may be empty (the origin), a scalar, one row, or a row per row of
// X. `sigma` may be empty (the identity), a d x d covariance, or a 1 x d
// ROW, which MATLAB reads as the DIAGONAL of a covariance and not as a
// 1 x d covariance -- the shape alone decides, so a one-dimensional call
// with two points cannot spell its variance as a row.
static bool normalDensity(const ICoreMatrix& x, const ICoreMatrix& mu,
const ICoreMatrix& sigma, const bool& haveMu,
const bool& haveSigma, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `mvtpdf(X, C, df)` -- the multivariate Student t density. `C` is a
// correlation matrix; see the header note on what MATLAB does when it is
// handed a covariance instead.
static bool studentDensity(const ICoreMatrix& x, const ICoreMatrix& correlation,
const ICoreMatrix& degrees, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `mvncdf(XU[, mu, Sigma])` and `mvncdf(XL, XU, mu, Sigma)` -- the
// probability in the box, one per row, as a column.
//
// `lower` empty is MATLAB's three-argument form, whose lower limit is
// -Inf in every coordinate. The refusal above three columns is this row's
// and names Z11.
static bool normalCumulative(const ICoreMatrix& lower, const ICoreMatrix& upper,
const ICoreMatrix& mu, const ICoreMatrix& sigma,
const bool& haveLower, const bool& haveMu,
const bool& haveSigma, ICoreMatrix& out,
std::string* whyNot = nullptr);
// The bivariate standardized cdf itself, exposed because two rows want it
// and one implementation is the point: `internal.stats.bvncdf(b, rho)`,
// P(Z1 <= b1, Z2 <= b2) for the standard bivariate normal of correlation
// rho. `b` entries may be -Inf.
static double bivariateNormalCdf(const double& first, const double& second,
const double& rho);
// ---- the four samplers ------------------------------------------------
//
// ⚠ None of these four can be parity-compared, for the reason T3.3's
// `normrnd` states: MATLAB's `randn` is a Mersenne Twister stream and this
// console's is not, so the same seed is not the same draw and no tolerance
// band makes it one. They cross with `Compare::Skip` (C0.9, board rule
// Z11) and are pinned by regress cases on a derived property -- a shape, a
// count, a symmetry -- exactly as `rand`'s own cases are.
// `mvnrnd(mu, Sigma[, n])` -- n rows, each `mu + z * chol(Sigma)`.
static bool normalRandom(const ICoreMatrix& mu, const ICoreMatrix& sigma,
const size_t& rows, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `mvtrnd(C, df[, n])` -- n rows of the multivariate t, each a normal row
// divided by the square root of its own scaled chi-square draw.
static bool studentRandom(const ICoreMatrix& correlation, const double& degrees,
const size_t& rows, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `wishrnd(Sigma, df)` and `iwishrnd(Sigma, df)` -- one d x d draw.
// Both go through Smith and Hocking's Bartlett decomposition, which is
// what `wishrnd.m` does: a lower triangle of standard normals under a
// diagonal of chi roots, so a draw costs d(d+1)/2 numbers rather than a
// df x d sample.
static bool wishartRandom(const ICoreMatrix& sigma, const double& degrees,
const bool& inverse, ICoreMatrix& out,
std::string* whyNot = nullptr);
};
};ICorePenalizedRegression.h#
src/ICoreBlocks/ICoreMath/Statistics/ICorePenalizedRegression.h
ICorePenalizedRegression#
ICorePenalizedRegression.h:48 · class · nested LassoOptions · 0 declaration(s)
The penalized and latent-variable regressions of toolbox row T3.15: lasso, ridge and plsregress.
class ICorePenalizedRegression {
public:
// `lasso`'s name/value options, with `lasso.m`'s own defaults. The four
// that are NOT here are refused at the console rather than defaulted:
// `'CV'` and `'MCReps'` (random folds -- Z11 -- and a `cvpartition` is an
// object, Z3), `'Options'` (a parallel-pool struct) and
// `'PredictorNames'` (a cellstr that only ever reaches the refused
// `FitInfo`).
struct LassoOptions {
// The elastic-net mixing, `0 < alpha <= 1`. 1 is the lasso; below it
// the ridge half of the penalty enters through `shrinkFactor`.
double alpha = 1.0;
// The default sequence: `numLambda` values geometric from `lambdaMax`
// down to `lambdaMax * lambdaRatio`. ⚠ `lambdaRatio == 0` is not "no
// ratio" -- it means "use 1e-4 and then make the LAST value exactly
// zero", which is the unpenalized fit and is how `lasso.m` spells it.
size_t numLambda = 100;
double lambdaRatio = 1.0e-4;
// A caller-supplied sequence, sorted DESCENDING before use. When this
// is non-empty the overfit break below is disabled, because a caller
// who named the penalties gets all of them.
std::vector<double> lambda;
// `'DFmax'`: stop once a fit has more than this many nonzero
// coefficients. 0 means "no limit" (`lasso.m` clamps it to `P`).
size_t degreesOfFreedomMax = 0;
bool standardize = true;
double relativeTolerance = 1.0e-4;
size_t maximumIterations = 100000;
bool intercept = true;
// `'Weights'`, one per observation; normalised to sum to 1 as
// `standardizeW` does. Empty is unweighted, which is NOT the same
// arithmetic as equal weights -- the unweighted path divides by `N`
// where the weighted one divides by 1.
std::vector<double> weights;
// `'UseCovariance'`: -1 is MATLAB's `'auto'` (false when `N <= D`,
// otherwise true unless the Gram matrix would exceed `cacheSize`
// megabytes). The two paths are the same algebra reached two ways --
// one keeps the residual, the other the Gram matrix -- so they differ
// in the last digits and BOTH sides must take the same branch.
int useCovariance = -1;
double cacheSize = 1.0e3;
};
// `B = lasso(X, y, ...)`. `coefficients` is `P x nLambda` on the ORIGINAL
// scale of `X`, its columns in ASCENDING lambda order -- `lasso.m` fits
// from `lambdaMax` downward and reverses everything at the end, so column
// 1 is the least penalized fit.
//
// The three other outputs are the numeric fields of the `FitInfo` struct
// this console has no kind for (T0.6): `lambdaUsed` (1 x nLambda),
// `intercepts` (1 x nLambda) and `meanSquaredError` (1 x nLambda). They
// are filled here so that the row lands whole and the console can answer
// them the day a record kind exists; nothing but the coefficients crosses
// the bridge today.
//
// ⚠ The sequence can come back SHORTER than `numLambda`: `lasso.m` stops
// early when a fit exceeds `DFmax`, and -- only for a sequence it built
// itself -- when the mean squared error falls below `1e-3` of the null
// model's, which it reads as overfitting.
static bool lasso(const ICoreMatrix& x, const ICoreMatrix& y, const LassoOptions& options,
ICoreMatrix& coefficients, ICoreMatrix& lambdaUsed,
ICoreMatrix& intercepts, ICoreMatrix& meanSquaredError,
std::string* whyNot = nullptr);
// `b = ridge(y, X, k[, scaled])`. `k` may hold several penalties, one per
// COLUMN of the answer. `restoreScale` is MATLAB's `scaled == 0`: see the
// warning at the top of this file -- the default is the scaled fit.
//
// MATLAB solves it as one tall least squares, `[Z; sqrt(k) I] \ [y; 0]`,
// rather than by forming `(Z'Z + kI)`, and so does this: the two are the
// same matrix in exact arithmetic and the augmented one is the better
// conditioned of the two by a square root.
static bool ridge(const ICoreMatrix& y, const ICoreMatrix& x, const ICoreMatrix& penalties,
const bool& restoreScale, ICoreMatrix& coefficients,
std::string* whyNot = nullptr);
// `[XL, YL, XS, YS, beta, pctVar, MSE] = plsregress(X, Y[, ncomp])`,
// SIMPLS exactly as `plsregress.m`'s local `simpls` writes it, including
// the TWICE-repeated Gram-Schmidt of the deflation basis and of the Y
// scores -- the repetition is not belt and braces, it is what keeps the
// basis orthogonal once the covariance has been deflated a few times.
//
// `components` of 0 means `min(n - 1, dx)`, MATLAB's default. `beta` is
// `(dx + 1) x dy` with the intercept row when `intercept` is true and
// `dx x dy` without it. The eighth output `stats` is a struct and is
// refused by name (T0.6).
static bool partialLeastSquares(const ICoreMatrix& x, const ICoreMatrix& y,
const size_t& components, const bool& intercept,
ICoreMatrix& xLoadings, ICoreMatrix& yLoadings,
ICoreMatrix& xScores, ICoreMatrix& yScores,
ICoreMatrix& beta, ICoreMatrix& percentVariance,
ICoreMatrix& meanSquaredError, std::string* whyNot = nullptr);
};
};ICorePrincipalComponents.h#
src/ICoreBlocks/ICoreMath/Statistics/ICorePrincipalComponents.h
ICorePrincipalComponents#
ICorePrincipalComponents.h:41 · class · 0 declaration(s)
pca, pcacov and pcares -- toolbox row T3.17.
class ICorePrincipalComponents {
public:
// `[coeff, score, latent, tsquared, explained, mu] = pca(X)`. `centred`
// is MATLAB's 'Centered' (default true); `components` is 'NumComponents'
// (0 for all), and dropping components changes only how many COLUMNS of
// coeff and score come back -- never latent or explained.
static bool principalComponents(const ICoreMatrix& x, const bool& centred,
const size_t& components, ICoreMatrix& coefficients,
ICoreMatrix& scores, ICoreMatrix& latent,
ICoreMatrix& tSquared, ICoreMatrix& explained,
ICoreMatrix& means, std::string* whyNot = nullptr);
// `[coeff, latent, explained] = pcacov(C)` -- the same decomposition of a
// COVARIANCE matrix instead of the data, through its SVD, with the same
// sign convention.
static bool principalComponentsOfCovariance(const ICoreMatrix& c, ICoreMatrix& coefficients,
ICoreMatrix& latent, ICoreMatrix& explained,
std::string* whyNot = nullptr);
// `[residuals, reconstructed] = pcares(X, ndim)` -- what the first `ndim`
// components do NOT explain. `ndim` above the usable rank is clamped to
// `min(n-1, p)`, as pcares.m clamps it, rather than refused.
static bool principalComponentResiduals(const ICoreMatrix& x, const size_t& dimensions,
ICoreMatrix& residuals, ICoreMatrix& reconstructed,
std::string* whyNot = nullptr);
};
};ICoreRegressionFits.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreRegressionFits.h
ICoreRegressionFits#
ICoreRegressionFits.h:34 · class · 1 declaration(s)
The regression fits of toolbox row T3.14: robustfit, glmfit, glmval and nlparci.
class ICoreRegressionFits {
public:
// `robustfit`'s nine weight functions, with the tuning constants
// `statrobustwfun.m` pairs them with. The constant is not a default in the
// ordinary sense -- each is the value that makes that weight function 95%
// as efficient as least squares on normal data, so changing one changes
// what the estimator IS.
enum class WeightFunction {
Andrews, Bisquare, Cauchy, Fair, Huber, Logistic, Ols, Talwar, Welsch,
};
static double defaultTuning(const WeightFunction& weight);
// MATLAB's five GLM distributions and the eight links. `Canonical` means
// "whichever link this distribution's row of `getGLMVariance.m` names",
// which is identity for the normal, logit for the binomial, log for the
// Poisson, reciprocal for the gamma and the -2 power for the inverse
// Gaussian -- so the canonical link is a property of the DISTRIBUTION and
// not a link of its own.
enum class GlmDistribution { Normal, Binomial, Poisson, Gamma, InverseGaussian };
enum class Link {
Canonical, Identity, Log, Logit, Probit, ComplementaryLogLog, LogLog, Reciprocal, Power,
};
// `b = robustfit(X, y[, wfun, tune, const])`. A constant column is ADDED
// unless `addConstant` is false -- the opposite of `regress`, which never
// adds one, and the trap the two functions share a page for.
static bool robustFit(const ICoreMatrix& x, const ICoreMatrix& y,
const WeightFunction& weight, const double& tune,
const bool& addConstant, ICoreMatrix& coefficients,
std::string* whyNot = nullptr);
// `b = glmfit(X, y, distr[, "link", L])`. `trials` is the binomial's N; it
// is ignored for every other distribution, as MATLAB ignores it with a
// warning. `linkPower` is the exponent when `link` is `Power`.
static bool generalizedLinearFit(const ICoreMatrix& x, const ICoreMatrix& y,
const GlmDistribution& distribution, const Link& link,
const double& linkPower, const ICoreMatrix& trials,
const bool& addConstant, ICoreMatrix& coefficients,
std::string* whyNot = nullptr);
// `yhat = glmval(b, X, link[, "size", N])` -- the inverse link of `X*b`,
// times N for a binomial. No fitting happens here; this is the forward
// half on its own, which is why it takes a link name rather than a
// distribution.
static bool generalizedLinearValue(const ICoreMatrix& coefficients, const ICoreMatrix& x,
const Link& link, const double& linkPower,
const ICoreMatrix& trials, const bool& addConstant,
ICoreMatrix& fitted, std::string* whyNot = nullptr);
// The model `nlinfit` fits, as the evaluator hands it over: the
// coefficient vector and the predictor matrix in, the fitted values out,
// and a reason when the handle cannot be called. A console function
// handle lives in the expression evaluator and not in this zone, which is
// the same seam `bootstrp`'s statistic crosses.
using Model = std::function<bool(const ICoreMatrix& beta, const ICoreMatrix& x,
ICoreMatrix& fitted, std::string& whyNot)>;
// `[beta, r, J, CovB, MSE] = nlinfit(X, y, modelfun, beta0)`.
//
// ⚠ This one IS the situation T3.6 met with `fzero`, and its note says so
// rather than hoping otherwise: Levenberg-Marquardt stops on `TolX` 1e-8
// of the step or `TolFun` 1e-8 of the relative change in the sum of
// squares, whichever comes first, and MATLAB's answer is where that
// stopped. So this transcribes the ITERATION -- the 0.01 starting damping,
// the tenfold increase on a failed step and the tenth on a good one, the
// augmented `[J; sqrt(lambda*diag(J'J))]` least-squares solve, and
// `statjacobian`'s forward difference at `eps^(1/3)` times each
// coefficient -- rather than solving the problem its own way.
static bool nonlinearFit(const ICoreMatrix& x, const ICoreMatrix& y, const Model& model,
const ICoreMatrix& start, ICoreMatrix& coefficients,
ICoreMatrix& residual, ICoreMatrix& jacobian,
ICoreMatrix& covariance, double& meanSquaredError,
std::string* whyNot = nullptr);
// `[ypred, delta] = nlpredci(modelfun, X, beta, resid, ...)`. `newObservation`
// is MATLAB's `"predopt", "observation"` (the interval for a future
// OBSERVATION, which carries the error variance as well as the parameter
// uncertainty); `simultaneous` is its `"simopt", "on"`.
//
// ⚠ Its finite-difference step is NOT `statjacobian`'s. Where `nlinfit`
// falls back to `DerivStep * norm(beta)` for a coefficient that is zero,
// `nlpredci.m` falls back to `DerivStep * sqrt(norm(beta))` -- the square
// root of the norm. Two files, two rules, and the difference only shows
// on a model fitted from a zero start.
static bool nonlinearPredictionInterval(const Model& model, const ICoreMatrix& x,
const ICoreMatrix& beta, const ICoreMatrix& residual,
const ICoreMatrix& matrix, const bool& isCovariance,
const double& alpha, const bool& newObservation,
const bool& simultaneous, ICoreMatrix& predicted,
ICoreMatrix& halfWidth, std::string* whyNot = nullptr);
// `ci = nlparci(beta, resid, "covar", C)` or `("jacobian", J)`. Two ways
// in, one answer: a standard error per coefficient, times `tinv` on
// `n - p` degrees of freedom. `isCovariance` says which of the two the
// matrix is.
static bool nonlinearParameterInterval(const ICoreMatrix& beta, const ICoreMatrix& residual,
const ICoreMatrix& matrix, const bool& isCovariance,
const double& alpha, ICoreMatrix& interval,
std::string* whyNot = nullptr);
// `[b, se, pval, inmodel] = stepwisefit(X, y[, "penter", a][, "premove", b])`.
//
// ⚠ The answer is an ORDER of decisions, not a formula: at each step the
// out-of-model column with the smallest p-value enters if that p is below
// `penter`, otherwise the in-model column with the largest p leaves if
// that p is above `premove`, and the loop stops when neither applies. So
// what has to agree with MATLAB is the SEQUENCE, and the two p-values it
// compares are computed on different degrees of freedom -- `dfe` for a
// column in the model and `dfe - 1` for one outside it, because entering
// it would spend one more.
//
// `scaled` is MATLAB's `"scale"`: the columns are standardised for the
// arithmetic either way, and this says whether the coefficients are
// reported on that scale or divided back.
static bool stepwiseFit(const ICoreMatrix& x, const ICoreMatrix& y, const double& enter,
const double& remove, const bool& scaled,
ICoreMatrix& coefficients, ICoreMatrix& standardError,
ICoreMatrix& probability, ICoreMatrix& inModel,
std::string* whyNot = nullptr);
};
};ICoreStatisticalResampling.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreStatisticalResampling.h
ICoreStatisticalResampling#
ICoreStatisticalResampling.h:31 · class · 1 declaration(s)
Toolbox row T3.23 -- jackknife, bootstrp, bootci, randsample and datasample: the five ways the Statistics Toolbox draws a new sample out of an old one.
class ICoreStatisticalResampling {
public:
// The statistic, as the evaluator hands it over: the resampled inputs in,
// the statistic flattened to a row out, and a reason when it cannot.
using Statistic = std::function<bool(const std::vector<ICoreMatrix>&,
std::vector<double>&, std::string&)>;
// `jackstat = jackknife(jackfun, X1, ...)` -- one ROW per left-out
// observation, in order, so row i is the statistic of the sample without
// row i. Deterministic, and the only name here that is.
static bool leaveOneOut(const std::vector<ICoreMatrix>& data, const Statistic& statistic,
ICoreMatrix& out, std::string* whyNot = nullptr);
// `bootstat = bootstrp(nboot, bootfun, d1, ...)` -- `nboot` rows, each the
// statistic of a sample of n rows drawn WITH replacement.
static bool bootstrap(const size_t& draws, const std::vector<ICoreMatrix>& data,
const Statistic& statistic, ICoreMatrix& out,
std::string* whyNot = nullptr);
// MATLAB's `bootci` interval types. `"student"` is refused by name: it is a
// bootstrap inside a bootstrap, which is a different amount of computation
// and a different row.
enum class Interval { BiasCorrectedAccelerated, Percentile, Normal };
static bool intervalOf(const std::string& word, Interval& type);
// `ci = bootci(nboot, bootfun, d1, ...)` -- a 2 x k interval, the lower row
// then the upper. The default type is the bias-corrected and accelerated
// one, whose acceleration term is a JACKKNIFE of the same statistic, so a
// `bootci` costs a `bootstrp` plus a `jackknife`.
static bool bootstrapInterval(const size_t& draws, const std::vector<ICoreMatrix>& data,
const Statistic& statistic, const double& alpha,
const Interval& type, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `y = randsample(n, k[, replace[, w]])` and `randsample(population, k, ...)`.
// `population` empty means the first form, where the population is 1..n.
static bool randomSample(const double& count, const ICoreMatrix& population,
const size_t& draws, const bool& withReplacement,
const ICoreMatrix& weights, ICoreMatrix& out,
std::string* whyNot = nullptr);
// `y = datasample(data, k[, dim][, "Replace", tf][, "Weights", w])` -- the
// same draw as `randsample`, applied to the ROWS (dim 1) or the COLUMNS
// (dim 2) of a matrix. The default is WITH replacement, which is the
// opposite of `randsample`'s default and is MATLAB's own inconsistency.
static bool dataSample(const ICoreMatrix& data, const size_t& draws,
const size_t& dimension, const bool& withReplacement,
const ICoreMatrix& weights, ICoreMatrix& out,
std::string* whyNot = nullptr);
};
};ICoreStatisticsPlots.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreStatisticsPlots.h
ICoreStatisticsPlots#
ICoreStatisticsPlots.h:32 · class · nested ProbabilityPlot, BoxSummary, HistogramFit, StepCurve · 3 declaration(s)
The NUMBERS behind the Statistics and Machine Learning Toolbox's plots (toolbox board row T3.25): histfit, normplot, probplot, qqplot, cdfplot and boxplot.
class ICoreStatisticsPlots {
public:
// `normplot` / `probplot` / `qqplot` against the normal: the sorted sample
// against the normal quantiles of its plotting positions, plus the
// reference line MATLAB draws through the first and third quartiles of
// both. `lineX`/`lineY` are the two ENDS of that line, extended to the
// data range the way `normplot.m` extends it.
struct ProbabilityPlot {
std::vector<double> sorted; // the sample, ascending
std::vector<double> quantiles; // norminv((i - 0.5)/n)
double lineFirstX = 0.0, lineFirstY = 0.0; // the reference line's ends
double lineLastX = 0.0, lineLastY = 0.0;
};
static bool probabilityPlot(const ICoreMatrix& x, ProbabilityPlot& out, std::string& whyNot);
// `qqplot(x, y)`: two samples against each other. The SHORTER sample sets
// the positions and the longer one is read at those percentiles, which is
// what `qqplot.m` does with `prctile` -- so two samples of different
// length still answer one curve.
static bool quantileQuantile(const ICoreMatrix& x, const ICoreMatrix& y,
ProbabilityPlot& out, std::string& whyNot);
// `boxplot`'s five numbers for one group, its notch and its outliers.
// ⚠ `whiskerLow`/`whiskerHigh` are the ADJACENT VALUES (observations),
// and they are clamped INTO the box when the most extreme non-outlier
// falls between the quartiles -- `boxplot.m`'s own `wlo = min(wlo, p25)`,
// which happens on small samples and looks like a bug when it is not.
struct BoxSummary {
double q1 = 0.0, median = 0.0, q3 = 0.0;
double whiskerLow = 0.0, whiskerHigh = 0.0;
double notchLow = 0.0, notchHigh = 0.0;
std::vector<double> outliers;
};
static bool boxSummary(const ICoreMatrix& x, BoxSummary& out, std::string& whyNot);
// `histfit`: the histogram MATLAB would draw (its own bin picker, at
// `ceil(sqrt(n))` bins unless told otherwise) and the fitted density,
// scaled to the histogram's AREA -- `n * binWidth * pdf`, which is what
// makes the curve sit on the bars rather than under them.
//
// The curve is 100 points over the fitted distribution's three-sigma
// range, `icdf(pd, [0.0013499 0.99865])` -- MATLAB's own two constants.
struct HistogramFit {
std::vector<double> edges;
std::vector<double> counts;
std::vector<double> curveX;
std::vector<double> curveY;
double mean = 0.0;
double deviation = 0.0;
};
static bool histogramFit(const ICoreMatrix& x, const size_t& bins, HistogramFit& out,
std::string& whyNot);
// `cdfplot`: the empirical CDF as the STAIRS MATLAB draws -- each distinct
// value twice, so the curve rises vertically at the observation and runs
// flat between two. MATLAB pads the ends with -Inf and +Inf; this answers
// the finite part and the caller draws the flat ends it wants.
struct StepCurve {
std::vector<double> x;
std::vector<double> y;
};
static bool empiricalCdfCurve(const ICoreMatrix& x, StepCurve& out, std::string& whyNot);
};
};ICoreSummaryStatistics.h#
src/ICoreBlocks/ICoreMath/Statistics/ICoreSummaryStatistics.h
ICoreSummaryStatistics#
ICoreSummaryStatistics.h:52 · class · 1 declaration(s)
The six summary statistics of toolbox row T3.2 -- geomean, harmmean, trimmean, tabulate, crosstab and grpstats.
class ICoreSummaryStatistics {
public:
// `geomean(x[, dim])` -- exp(mean(log(x))). A NEGATIVE entry is an error
// on both sides (geomean.m raises stats:geomean:BadData); a zero is not,
// and answers 0.
static ICoreMatrix geometricMean(const ICoreMatrix& x, const int& dim,
std::string* whyNot = nullptr);
// `harmmean(x[, dim])` -- 1 ./ mean(1 ./ x). harmmean.m checks nothing, so
// neither does this: a zero entry answers 0 (1/0 is Inf, its mean is Inf,
// and 1/Inf is 0) and a sign mixture answers whatever the arithmetic gives.
static ICoreMatrix harmonicMean(const ICoreMatrix& x, const int& dim,
std::string* whyNot = nullptr);
// `trimmean(x, percent[, flag[, dim]])`. The flag chooses how many samples
// come off each end; see the class comment for why they differ.
enum class Trim { Round, Floor, Weighted };
static ICoreMatrix trimmedMean(const ICoreMatrix& x, const double& percent,
const Trim& flag, const int& dim,
std::string* whyNot = nullptr);
// `tabulate(x)` -- the [value, count, percent] table of a numeric VECTOR,
// NaN dropped. Two answer shapes; see the class comment.
static ICoreMatrix tabulate(const ICoreMatrix& x, std::string* whyNot = nullptr);
// `[table, chi2, p] = crosstab(x, y)` -- the contingency table of two
// grouping vectors of the same length, with Pearson's statistic and its
// upper-tail p-value. `chi2` and `p` are NaN when either variable has a
// single level, which is MATLAB's answer for a table that cannot be
// independent or dependent.
static bool crossTabulation(const ICoreMatrix& x, const ICoreMatrix& y,
ICoreMatrix& table, double& chiSquare, double& pValue,
std::string* whyNot = nullptr);
// `[means, sems, counts] = grpstats(x, group)` -- x's rows grouped by the
// corresponding entry of `group`, one output row per group in ascending
// group order. `sems` is `std(x, 0, 1) / sqrt(n)` per group, so a group of
// one answers 0 and not NaN (MATLAB's std of one sample is 0).
//
// `counts` has x's WIDTH rather than one column -- a 4 x 2 x over two
// groups answers a 2 x 2 of group sizes, not a 2 x 1 -- because MATLAB
// applies its count function per column. Measured, not assumed.
//
// MATLAB's fourth output is `gname`, a CELL of group labels; the console
// has no cell kind, so the caller is told to read the group values out of
// `unique(group)` instead (T0.6, Z3).
static bool groupStatistics(const ICoreMatrix& x, const ICoreMatrix& group,
ICoreMatrix& means, ICoreMatrix& sems, ICoreMatrix& counts,
std::string* whyNot = nullptr);
};
};