Report a bug
If you spot a problem with this page, click here to create a GitHub issue.
Improve this page
Quickly fork, edit online, and submit a pull request for this page. Requires a signed-in GitHub account. This works well for small changes. If you'd like to make larger changes you may want to consider using a local clone.

mir.stat.descriptive

This module contains algorithms for descriptive statistics.
Authors:
John Michael Hall, Ilya Yaroshenko
enum QuantileAlgo: int;
Algorithms used to calculate the quantile of an input x at probability p.
These algorithms match the same provided in R's (as of version 3.6.2) quantile function. In turn, these were discussed in Hyndman and Fan (1996).
All sample quantiles are defined as weighted averages of consecutive order statistics. For each QuantileAlgo, the sample quantile is given by (using R's 1-based indexing notation):
(1 - gamma) * x_{j} + gamma * x_{j + 1}
where x_{j} is the jth order statistic. gamma is a function of j = floor(np + m) and g = np + m - j where n is the sample size, p is the probability, and m is a constant determined by the quantile type.

Discontinuous sample quantile

Type m gamma
type1 0, 0 if g = 0 and 1 otherwise.
type2 0, 0.5 if g = 0 and 1 otherwise.
type3 -0.5, 0 if g = 0 and j is even and 1 otherwise.

Continuous sample quantile

, Type m gamma type4 0, gamma = g type5 0.5, gamma = g type6 p, gamma = g type7 1 - p, gamma = g type8 (p + 1) / 3, gamma = g type9 p / 4 + 3 / 8, gamma = g

References Hyndman, R. J. and Fan, Y. (1996) Sample quantiles in statistical packages, American Statistician 50, 361--365. 10.2307/2684934.

See Also:
type1

Discontinuous sample quantile

Inverse of empirical distribution function.
type2
Similar to type1, but averages at discontinuities.
type3
SAS definition: nearest even order statistic.
type4

Continuous sample quantile

Linear interpolation of the empirical cdf.
type5
A piece-wise linear function hwere the knots are the values midway through the steps of the empirical cdf. Popular amongst hydrologists.
type6
Used by Minitab and by SPSS.
type7
This is used by S and is the default for R.
type8
The resulting quantile estimates are approximately median-unbiased regardless of the distribution of the input. Preferred by Hyndman and Fan (1996).
type9
The resulting quantile estimates are approximately unbiased for the expected order statistics of the input is normally distributed.
template quantile(F, QuantileAlgo quantileAlgo = QuantileAlgo.type7, bool allowModifySlice = false, bool allowModifyProbability = false) if (isFloatingPoint!F || (quantileAlgo == QuantileAlgo.type1 || quantileAlgo == QuantileAlgo.type3))
Computes the quantile(s) of the input, given one or more probabilities p.
By default, if p is a Slice , built-in dynamic array, or type with asSlice, then the output type is a reference-counted copy of the input. A run-time parameter is provided to instead overwrite the input in-place.
For all QuantileAlgo except QuantileAlgo.type1 and QuantileAlgo.type3, by default, if F is not floating point type or complex type, then the result will have a double type if F is implicitly convertible to a floating point type or have a cdouble type if F is implicitly convertible to a complex type.
For QuantileAlgo.type1 and QuantileAlgo.type3, the return type is the elementType  of the input.
Parameters:
F controls type of output
quantileAlgo algorithm for calculating quantile (default: QuantileAlgo.type7)
allowModifySlice controls whether the input is modified in place, default is false
Returns:
The quantile of all the elements in the input at probability p.
See Also:
quantileType!(F, quantileAlgo) quantile(Iterator, size_t N, SliceKind kind, G)(Slice!(Iterator, N, kind) slice, G p)
if (isFloatingPoint!(Unqual!G));
Parameters:
Slice!(Iterator, N, kind) slice slice
G p probability
auto quantile(IteratorA, size_t N, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, N, kindA) slice, Slice!(IteratorB, 1, kindB) p)
if (isFloatingPoint!(elementType!(Slice!IteratorB)));

auto quantile(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice, scope const F[] p...)
if (isFloatingPoint!(elementType!(F[])));

quantileType!(F, quantileAlgo) quantile(G)(F[] array, G p)
if (isFloatingPoint!(Unqual!G));

auto quantile(G)(F[] array, G[] p)
if (isFloatingPoint!(Unqual!G));

auto quantile(T, G)(T withAsSlice, G p)
if (hasAsSlice!T && isFloatingPoint!(Unqual!G));

auto quantile(T, U)(T withAsSlice, U p)
if (hasAsSlice!T && hasAsSlice!U);
Parameters:
Slice!(IteratorA, N, kindA) slice slice
Slice!(IteratorB, 1, kindB) p probability slice
template quantile(QuantileAlgo quantileAlgo = QuantileAlgo.type7, bool allowModifySlice = false, bool allowModifyProbability = false)

template quantile(F, string quantileAlgo, bool allowModifySlice = false, bool allowModifyProbability = false)

template quantile(string quantileAlgo, bool allowModifySlice = false, bool allowModifyProbability = false)
quantileType!(Slice!Iterator, quantileAlgo) quantile(Iterator, size_t N, SliceKind kind, G)(Slice!(Iterator, N, kind) slice, G p)
if (isFloatingPoint!(Unqual!G));

auto quantile(IteratorA, size_t N, SliceKind kindA, IteratorB, SliceKind kindB)(Slice!(IteratorA, N, kindA) slice, Slice!(IteratorB, 1, kindB) p)
if (isFloatingPoint!(elementType!(Slice!IteratorB)));

auto quantile(Iterator, size_t N, SliceKind kind, G)(Slice!(Iterator, N, kind) slice, scope G[] p...)
if (isFloatingPoint!(elementType!(G[])));

auto quantile(T, G)(T[] array, G p)
if (isFloatingPoint!(Unqual!G));

auto quantile(T, G)(T[] array, G[] p)
if (isFloatingPoint!(Unqual!G));

auto quantile(T, G)(T withAsSlice, G p)
if (hasAsSlice!T && isFloatingPoint!(Unqual!G));

auto quantile(T, U)(T withAsSlice, U p)
if (hasAsSlice!T && hasAsSlice!U);
Parameters:
Slice!(Iterator, N, kind) slice slice
G p probability
template interquartileRange(F, QuantileAlgo quantileAlgo = QuantileAlgo.type7, bool allowModifySlice = false)

template interquartileRange(QuantileAlgo quantileAlgo = QuantileAlgo.type7, bool allowModifySlice = false)

template interquartileRange(F, string quantileAlgo, bool allowModifySlice = false)

template interquartileRange(string quantileAlgo, bool allowModifySlice = false)
Computes the interquartile range of the input.
This function computes the result using quantile, i.e. result = quantile(x, 0.75) - quantile(x, 0.25).
For all QuantileAlgo except QuantileAlgo.type1 and QuantileAlgo.type3, by default, if F is not floating point type or complex type, then the result will have a double type if F is implicitly convertible to a floating point type or have a cdouble type if F is implicitly convertible to a complex type.
For QuantileAlgo.type1 and QuantileAlgo.type3, the return type is the elementType  of the input.
Parameters:
F controls type of output
quantileAlgo algorithm for calculating quantile (default: QuantileAlgo.type7)
allowModifySlice controls whether the input is modified in place, default is false
Returns:
The interquartile range of the input.
See Also:
Examples:
Simple example
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [3.0, 1.0, 4.0, 2.0, 0.0].sliced;

assert(x.interquartileRange.approxEqual(2.0));
Examples:
Interquartile Range of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [1.0, 9.8, 0.2, 8.5, 5.8, 3.5, 4.5, 8.2, 5.2, 5.2,
          2.5, 1.8, 2.2, 3.8, 5.2, 9.2, 6.2, 9.2, 9.2, 8.5].sliced;

assert(x.interquartileRange.approxEqual(5.25));
Examples:
Interquartile Range of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.slice: sliced;

auto x = [
    [1.0, 9.8, 0.2, 8.5, 5.8, 3.5, 4.5, 8.2, 5.2, 5.2],
    [2.5, 1.8, 2.2, 3.8, 5.2, 9.2, 6.2, 9.2, 9.2, 8.5]
].fuse;

assert(x.interquartileRange.approxEqual(5.25));
Examples:
Allow modification of input
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [3.0, 1.0, 4.0, 2.0, 0.0].sliced;
auto x_copy = x.dup;

auto result = x.interquartileRange!(QuantileAlgo.type7, true);
assert(!x.all!approxEqual(x_copy));
Examples:
Can also set algorithm type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [1.0, 9.8, 0.2, 8.5, 5.8, 3.5, 4.5, 8.2, 5.2, 5.2,
          2.5, 1.8, 2.2, 3.8, 5.2, 9.2, 6.2, 9.2, 9.2, 8.5].sliced;

assert(x.interquartileRange!"type1".approxEqual(6.0));
assert(x.interquartileRange!"type2".approxEqual(5.5));
assert(x.interquartileRange!"type3".approxEqual(6.0));
assert(x.interquartileRange!"type4".approxEqual(6.0));
assert(x.interquartileRange!"type5".approxEqual(5.5));
assert(x.interquartileRange!"type6".approxEqual(5.75));
assert(x.interquartileRange!"type7".approxEqual(5.25));
assert(x.interquartileRange!"type8".approxEqual(5.583333));
assert(x.interquartileRange!"type9".approxEqual(5.5625));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto a = [1, 1e34, 1, -1e34, 0].sliced;

auto x = a * 10_000;

auto result0 = x.interquartileRange!float;
assert(result0.approxEqual(10_000));
static assert(is(typeof(result0) == float));

auto result1 = x.interquartileRange!(float, "type8");
assert(result1.approxEqual(6.666667e37));
static assert(is(typeof(result1) == float));
Examples:
Support for array
import mir.math.common: approxEqual;

double[] x = [3.0, 1.0, 4.0, 2.0, 0.0];
          
assert(x.interquartileRange.approxEqual(2.0));
quantileType!(F, quantileAlgo) interquartileRange(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);
Parameters:
Slice!(Iterator, N, kind) slice slice
quantileType!(F[], quantileAlgo) interquartileRange(scope F[] array...);
Parameters:
F[] array array
auto interquartileRange(T)(T withAsSlice)
if (hasAsSlice!T);
Parameters:
T withAsSlice withAsSlice
template medianAbsoluteDeviation(F)

meanType!(Slice!(Iterator, N, kind)) medianAbsoluteDeviation(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);

meanType!(T[]) medianAbsoluteDeviation(T)(scope const T[] ar...);

auto medianAbsoluteDeviation(T)(T withAsSlice)
if (hasAsSlice!T);
Calculates the median absolute deviation about the median of the input.
By default, if F is not floating point type, then the result will have a double type if F is implicitly convertible to a floating point type.
Parameters:
F output type
Returns:
The median absolute deviation of the input
Examples:
Simple example
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.medianAbsoluteDeviation.approxEqual(1.25));
Examples:
Median Absolute Deviation of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.medianAbsoluteDeviation.approxEqual(1.25));
Examples:
Median Absolute Deviation of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.medianAbsoluteDeviation.approxEqual(1.25));
Examples:
Median Absolute Deviation of dynamic array
import mir.math.common: approxEqual;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0];

assert(x.medianAbsoluteDeviation.approxEqual(1.25));
meanType!F medianAbsoluteDeviation(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);
Parameters:
Slice!(Iterator, N, kind) slice slice
template dispersion(alias centralTendency = mean, alias transform = "a * a", alias summarize = mean)
Calculates the dispersion of the input. For an input x, this function first centers x by subtracting each e in x by the result of centralTendency, then it transforms the centered values using the function transform, and then finally summarizes that information using the summarize funcion. The default functions provided are equivalent to calculating the population variance. The centralTendency default is the mean function, which results in the input being centered about the mean. The default transform function will square the centered values. The default summarize function is mean, which will return the mean of the squared centered values.
Parameters:
centralTendency function that will produce the value that the input is centered about, default is mean
transform function to transform centered values, default squares the centered values
summarize function to summarize the transformed centered values, default is mean
Returns:
The dispersion of the input
Examples:
Simple examples
import mir.math.common: approxEqual;
import mir.functional: naryFun;
import mir.ndslice.slice: sliced;

assert(dispersion([1.0, 2, 3]).approxEqual(2.0 / 3));

assert(dispersion([1.0 + 3i, 2, 3]).approxEqual((-4.0 - 6i) / 3));

assert(dispersion!(mean!float, "a * a", mean!float)([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(17.5 / 6));

static assert(is(typeof(dispersion!(mean!float, "a ^^ 2", mean!float)([1, 2, 3])) == float));
Examples:
Dispersion of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.dispersion.approxEqual(54.76562 / 12));
Examples:
Dispersion of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.dispersion.approxEqual(54.76562 / 12));
Examples:
Column dispersion of matrix
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: alongDim, byDim, map;

auto x = [
    [0.0,  1.0, 1.5, 2.0], 
    [3.5, 4.25, 2.0, 7.5],
    [5.0,  1.0, 1.5, 0.0]
].fuse;
auto result = [13.16667 / 3, 7.041667 / 3, 0.1666667 / 3, 30.16667 / 3];

// Use byDim or alongDim with map to compute dispersion of row/column.
assert(x.byDim!1.map!dispersion.all!approxEqual(result));
assert(x.alongDim!0.map!dispersion.all!approxEqual(result));

// FIXME
// Without using map, computes the dispersion of the whole slice
// assert(x.byDim!1.dispersion == x.sliced.dispersion);
// assert(x.alongDim!0.dispersion == x.sliced.dispersion);
Examples:
Can also set functions to change type of dispersion that is used
import mir.functional: naryFun;
import mir.math.common: approxEqual, fabs, sqrt;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;
          
alias square = naryFun!"a * a";

// Other population variance examples
assert(x.dispersion.approxEqual(54.76562 / 12));
assert(x.dispersion!mean.approxEqual(54.76562 / 12));
assert(x.dispersion!(mean, square).approxEqual(54.76562 / 12));
assert(x.dispersion!(mean, square, mean).approxEqual(54.76562 / 12));

// Population standard deviation
assert(x.dispersion!(mean, square, mean).sqrt.approxEqual(sqrt(54.76562 / 12)));

// Mean absolute deviation about the mean
assert(x.dispersion!(mean, fabs, mean).approxEqual(21.0 / 12));
//Mean absolute deviation about the median
assert(x.dispersion!(median, fabs, mean).approxEqual(19.25000 / 12));
//Median absolute deviation about the mean
assert(x.dispersion!(mean, fabs, median).approxEqual(1.43750));
//Median absolute deviation about the median
assert(x.dispersion!(median, fabs, median).approxEqual(1.25000));
Examples:
For integral slices, pass output type to centralTendency, transform, and summary functions as template parameter to ensure output type is correct.
import mir.functional: naryFun;
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

alias square = naryFun!"a * a";

auto y = x.dispersion;
assert(y.approxEqual(50.91667 / 12));
static assert(is(typeof(y) == double));

assert(x.dispersion!(mean!float, square, mean!float).approxEqual(50.91667 / 12));
Examples:
Dispersion works for complex numbers and other user-defined types (provided that the centralTendency, transform, and summary functions are defined for those types)
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [1.0 + 2i, 2 + 3i, 3 + 4i, 4 + 5i].sliced;
assert(x.dispersion.approxEqual((0.0+10.0i)/ 4));
Examples:
Compute mean tensors along specified dimention of tensors
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: as, iota, alongDim, map, repeat;

auto x = [
    [0.0, 1, 2],
    [3.0, 4, 5]
].fuse;

assert(x.dispersion.approxEqual(17.5 / 6));

auto m0 = [2.25, 2.25, 2.25];
assert(x.alongDim!0.map!dispersion.all!approxEqual(m0));
assert(x.alongDim!(-2).map!dispersion.all!approxEqual(m0));

auto m1 = [2.0 / 3, 2.0 / 3];
assert(x.alongDim!1.map!dispersion.all!approxEqual(m1));
assert(x.alongDim!(-1).map!dispersion.all!approxEqual(m1));

assert(iota(2, 3, 4, 5).as!double.alongDim!0.map!dispersion.all!approxEqual(repeat(1800.0 / 2, 3, 4, 5)));
Examples:
Arbitrary dispersion
import mir.functional: naryFun;
import mir.math.common: approxEqual;

alias square = naryFun!"a * a";

assert(dispersion(1.0, 2, 3).approxEqual(2.0 / 3));
assert(dispersion!(mean!float, square, mean!float)(1, 2, 3).approxEqual(2f / 3));
auto dispersion(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);

auto dispersion(T)(scope const T[] ar...);

auto dispersion(T)(T withAsSlice)
if (hasAsSlice!T);
Parameters:
Slice!(Iterator, N, kind) slice slice
enum SkewnessAlgo: int;
Skew algorithms. See Also: Skewness. Algorithms for calculating variance.
online
Similar to Welford's algorithm for updating variance, but adjusted for skewness. Can also put another SkewnessAccumulator of the same type, which uses the parallel algorithm from Terriberry that extends the work of Chan et al.
naive
Calculates skewness using (E(x^^3) - 3 * mu * sigma ^^ 2 + mu ^^ 3) / (sigma ^^ 3) (alowing for adjustments for population/sample skewness). This algorithm can be numerically unstable.
twoPass
Calculates skewness using a two-pass algorithm whereby the input is first scaled by the mean and variance (using VarianceAccumulator.online ) and then the sum of cubes is calculated from that.
threePass
Calculates skewness using a three-pass algorithm whereby the input is first scaled by the mean and variance (using VarianceAccumulator.twoPass ) and then the sum of cubes is calculated from that.
assumeZeroMean
Calculates skewness assuming the mean of the input is zero.
struct SkewnessAccumulator(T, SkewnessAlgo skewnessAlgo, Summation summation) if (isMutable!T && (skewnessAlgo == SkewnessAlgo.naive));
Examples:
naive
import mir.math.common: approxEqual, pow;
import mir.math.sum: Summation;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;

SkewnessAccumulator!(double, SkewnessAlgo.naive, Summation.naive) v;
v.put(x);
assert(v.skewness(PopulationTrueRT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));

v.put(4.0);
assert(v.skewness(PopulationTrueRT).approxEqual((100.238166 / 13) / pow(57.019231 / 13, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((100.238166 / 13) / pow(57.019231 / 13, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((100.238166 / 13) / pow(57.019231 / 12, 1.5) * (13.0 ^^ 2) / (12.0 * 11.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((100.238166 / 13) / pow(57.019231 / 12, 1.5) * (13.0 ^^ 2) / (12.0 * 11.0)));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
VarianceAccumulator!(T, VarianceAlgo.naive, summation) varianceAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) sumOfCubes;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
@property F skewness(F = T)(bool isPopulation)
if (isFloatingPoint!F);
struct SkewnessAccumulator(T, SkewnessAlgo skewnessAlgo, Summation summation) if (isMutable!T && (skewnessAlgo == SkewnessAlgo.online));
Examples:
online
import mir.math.common: approxEqual, pow;
import mir.math.sum: Summation;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;

SkewnessAccumulator!(double, SkewnessAlgo.online, Summation.naive) v;
v.put(x);
assert(v.skewness(PopulationTrueRT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));

v.put(4.0);
assert(v.skewness(PopulationTrueRT).approxEqual((100.238166 / 13) / pow(57.019231 / 13, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((100.238166 / 13) / pow(57.019231 / 13, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((100.238166 / 13) / pow(57.019231 / 12, 1.5) * (13.0 ^^ 2) / (12.0 * 11.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((100.238166 / 13) / pow(57.019231 / 12, 1.5) * (13.0 ^^ 2) / (12.0 * 11.0)));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
MeanAccumulator!(T, summation) meanAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) centeredSumOfSquares;
Summator!(T, summation) centeredSumOfCubes;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
void put()(SkewnessAccumulator!(T, skewnessAlgo, summation) v);
@property F skewness(F = T)(bool isPopulation)
if (isFloatingPoint!F);
struct SkewnessAccumulator(T, SkewnessAlgo skewnessAlgo, Summation summation) if (isMutable!T && (skewnessAlgo == SkewnessAlgo.twoPass || skewnessAlgo == SkewnessAlgo.threePass));
Examples:
twoPass & threePass
import mir.math.common: approxEqual, sqrt;
import mir.math.sum: Summation;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;

auto v = SkewnessAccumulator!(double, SkewnessAlgo.twoPass, Summation.naive)(x);
assert(v.skewness(PopulationTrueRT).approxEqual(12.000999 / 12));
assert(v.skewness(PopulationTrueCT).approxEqual(12.000999 / 12));
assert(v.skewness(PopulationFalseRT).approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));
assert(v.skewness(PopulationFalseCT).approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));

auto w = SkewnessAccumulator!(double, SkewnessAlgo.threePass, Summation.naive)(x);
assert(w.skewness(PopulationTrueRT).approxEqual(12.000999 / 12));
assert(w.skewness(PopulationTrueCT).approxEqual(12.000999 / 12));
assert(w.skewness(PopulationFalseRT).approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));
assert(w.skewness(PopulationFalseCT).approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));
size_t count;
Summator!(T, summation) scaledSumOfCubes;
this(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);
this(U)(U[] array);
this(T)(T withAsSlice)
if (hasAsSlice!T);
@property F skewness(F = T)(bool isPopulation)
if (isFloatingPoint!F);
struct SkewnessAccumulator(T, SkewnessAlgo skewnessAlgo, Summation summation) if (isMutable!T && (skewnessAlgo == SkewnessAlgo.assumeZeroMean));
Examples:
assumeZeroMean
import mir.math.common: approxEqual, pow;
import mir.math.stat: center;
import mir.math.sum: Summation;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;
auto x = a.center;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;

SkewnessAccumulator!(double, SkewnessAlgo.assumeZeroMean, Summation.naive) v;
v.put(x);
assert(v.skewness(PopulationTrueRT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * 12.0 ^^ 2 / (11.0 * 10.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * 12.0 ^^ 2 / (11.0 * 10.0)));

v.put(4.0);
assert(v.skewness(PopulationTrueRT).approxEqual((181.005859 / 13) / pow(70.765625 / 13, 1.5)));
assert(v.skewness(PopulationTrueCT).approxEqual((181.005859 / 13) / pow(70.765625 / 13, 1.5)));
assert(v.skewness(PopulationFalseRT).approxEqual((181.005859 / 13) / pow(70.765625 / 12, 1.5) * 13.0 ^^ 2 / (12.0 * 11.0)));
assert(v.skewness(PopulationFalseCT).approxEqual((181.005859 / 13) / pow(70.765625 / 12, 1.5) * 13.0 ^^ 2 / (12.0 * 11.0)));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
VarianceAccumulator!(T, VarianceAlgo.assumeZeroMean, summation) varianceAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) centeredSumOfCubes;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
void put()(SkewnessAccumulator!(T, skewnessAlgo, summation) v);
@property F skewness(F = T)(bool isPopulation)
if (isFloatingPoint!F);
template skewness(F, SkewnessAlgo skewnessAlgo = SkewnessAlgo.online, Summation summation = Summation.appropriate)

template skewness(SkewnessAlgo skewnessAlgo = SkewnessAlgo.online, Summation summation = Summation.appropriate)

template skewness(F, string skewnessAlgo, string summation = "appropriate")

template skewness(string skewnessAlgo, string summation = "appropriate")
Calculates the skewness of the input
By default, if F is not floating point type, then the result will have a double type if F is implicitly convertible to a floating point type.
Parameters:
F controls type of output
skewnessAlgo algorithm for calculating skewness (default: SkewnessAlgo.online) summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The skewness of the input, must be floating point or complex type
Examples:
Simple example
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

assert(skewness([1.0, 2, 3]).approxEqual(0.0));

assert(skewness([1.0, 2, 4]).approxEqual((2.222222 / 3) / pow(4.666667 / 2, 1.5) * (3.0 ^^ 2) / (2.0 * 1.0)));
assert(skewness([1.0, 2, 4], true).approxEqual((2.222222 / 3) / pow(4.666667 / 3, 1.5)));

assert(skewness!float([0, 1, 2, 3, 4, 6].sliced(3, 2)).approxEqual(0.462910));

static assert(is(typeof(skewness!float([1, 2, 3])) == float));
Examples:
Skewness of vector
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.skewness.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));
Examples:
Skewness of matrix
import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.skewness.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));
Examples:
Column skewness of matrix
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: alongDim, byDim, map;

auto x = [
    [0.0,  1.0, 1.5, 2.0], 
    [3.5, 4.25, 2.0, 7.5],
    [5.0,  1.0, 1.5, 0.0]
].fuse;
auto result = [-1.090291, 1.732051, 1.732051, 1.229809];

// Use byDim or alongDim with map to compute skewness of row/column.
assert(x.byDim!1.map!skewness.all!approxEqual(result));
assert(x.alongDim!0.map!skewness.all!approxEqual(result));

// FIXME
// Without using map, computes the skewness of the whole slice
// assert(x.byDim!1.skewness == x.sliced.skewness);
// assert(x.alongDim!0.skewness == x.sliced.skewness);
Examples:
Can also set algorithm type
import mir.math.common: approxEqual, pow, sqrt;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto x = a + 100_000_000_000;

// The default online algorithm is numerically unstable in this case
auto y = x.skewness;
assert(!y.approxEqual((117.005859 / 12) / pow(54.765625 / 11, 1.5) * (12.0 ^^ 2) / (11.0 * 10.0)));

// The naive algorithm has an assert error in this case because standard
// deviation is calculated naively as zero. The skewness formula would then
// be dividing by zero. 
//auto z0 = x.skewness!(real, "naive");

// The two-pass algorithm is also numerically unstable in this case
auto z1 = x.skewness!"twoPass";
assert(!z1.approxEqual(12.000999 / 12 * sqrt(12.0 * 11.0) / 10.0));
assert(!z1.approxEqual(y));

// However, the three-pass algorithm is numerically stable in this case
auto z2 = x.skewness!"threePass";
assert(z2.approxEqual((12.000999 / 12) * sqrt(12.0 * 11.0) / 10.0));
assert(!z2.approxEqual(y));

// And the assumeZeroMean algorithm provides the incorrect answer, as expected
auto z3 = x.skewness!"assumeZeroMean";
assert(!z3.approxEqual(y));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

//Set population skewness, skewness algorithm, sum algorithm or output type

auto a = [1.0, 1e98, 1, -1e98].sliced;
auto x = a * 10_000;

bool populationTrueRT = true;
bool populationFalseRT = false;
enum PopulationTrueCT = true;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.skewness(populationFalseRT).approxEqual(0.0));
assert(x.skewness(populationTrueRT).approxEqual(0.0));
assert(x.skewness(PopulationTrueCT).approxEqual(0.0));

assert(x.skewness!("online").approxEqual(0.0));
assert(x.skewness!("online", "kbn").approxEqual(0.0));
assert(x.skewness!("online", "kb2").approxEqual(0.0));
assert(x.skewness!("online", "precise").approxEqual(0.0));
assert(x.skewness!(double, "online", "precise").approxEqual(0.0));
assert(x.skewness!(double, "online", "precise")(populationTrueRT).approxEqual(0.0));

auto y = [uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.skewness!(ulong, "threePass");
assert(z == 0.0);
static assert(is(typeof(z) == double));
Examples:
For integral slices, can pass output type as template parameter to ensure output type is correct.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.skewness;
assert(y.approxEqual(0.925493));
static assert(is(typeof(y) == double));

assert(x.skewness!float.approxEqual(0.925493));
Examples:
Skewness works for other user-defined types (provided they can be converted to a floating point)
import mir.ndslice.slice: sliced;

static struct Foo {
    float x;
    alias x this;
}

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f)];
assert(foo.skewness == 0f);
Examples:
Compute skewness along specified dimention of tensors
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: as, iota, alongDim, map, repeat;

auto x = [
    [0.0, 1, 3],
    [3.0, 4, 5],
    [6.0, 7, 7],
].fuse;

assert(x.skewness.approxEqual(-0.308571));

auto m0 = [0, 0.0, 0.0];
assert(x.alongDim!0.map!skewness.all!approxEqual(m0));
assert(x.alongDim!(-2).map!skewness.all!approxEqual(m0));

auto m1 = [0.935220, 0.0, -1.732051];
assert(x.alongDim!1.map!skewness.all!approxEqual(m1));
assert(x.alongDim!(-1).map!skewness.all!approxEqual(m1));
assert(iota(3, 4, 5, 6).as!double.alongDim!0.map!skewness.all!approxEqual(repeat(0.0, 4, 5, 6)));
Examples:
Arbitrary skewness
assert(skewness(1.0, 2, 3) == 0.0);
assert(skewness!float(1, 2, 3) == 0f);
stdevType!F skewness(Range)(Range r, bool isPopulation = false)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
bool isPopulation true if population skewness, false if sample skewness (default)
stdevType!F skewness(scope const F[] ar...);
Parameters:
F[] ar values
enum KurtosisAlgo: int;
Kurtosis algorithms. See Also: Kurtosis. Algorithms for calculating variance.
online
Similar to Welford's algorithm for updating variance, but adjusted for kurtosis. Can also put another KurtosisAccumulator of the same type, which uses the parallel algorithm from Terriberry that extends the work of Chan et al.
naive
Calculates kurtosis using (E(x^^4) - 4 * E(x) * E(x ^^ 3) + 6 * (E(x) ^^ 2) E(X ^^ 2) + 3 E(x) ^^ 4) / sigma ^ 2 (allowing for adjustments for population/sample kurtosis). This algorithm can be numerically unstable.
twoPass
Calculates kurtosis using a two-pass algorithm whereby the input is first scaled by the mean and variance (using VarianceAccumulator.online ) and then the sum of quarts is calculated from that.
threePass
Calculates kurtosis using a three-pass algorithm whereby the input is first scaled by the mean and variance (using VarianceAccumulator.twoPass ) and then the sum of quarts is calculated from that.
assumeZeroMean
Calculates kurtosis assuming the mean of the input is zero.
struct KurtosisAccumulator(T, KurtosisAlgo kurtosisAlgo, Summation summation) if (isMutable!T && (kurtosisAlgo == KurtosisAlgo.naive));
Examples:
naive
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;
enum RawTrueCT = true;
enum RawFalseCT = false;
bool RawTrueRT = true;
bool RawFalseRT = false;

KurtosisAccumulator!(double, KurtosisAlgo.naive, Summation.naive) v;
v.put(x);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));

v.put(4.0);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
MeanAccumulator!(T, summation) meanAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) sumOfSquares;
Summator!(T, summation) sumOfCubes;
Summator!(T, summation) sumOfQuarts;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
@property F kurtosis(F = T)(bool isPopulation, bool isRaw)
if (isFloatingPoint!F);
struct KurtosisAccumulator(T, KurtosisAlgo kurtosisAlgo, Summation summation) if (isMutable!T && (kurtosisAlgo == KurtosisAlgo.online));
Examples:
online
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;
enum RawTrueCT = true;
enum RawFalseCT = false;
bool RawTrueRT = true;
bool RawFalseRT = false;

KurtosisAccumulator!(double, KurtosisAlgo.online, Summation.naive) v;
v.put(x);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));

v.put(4.0);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual((745.608180 / 13) / pow(57.019231 / 13, 2.0) * (12.0 * 14.0) / (11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
MeanAccumulator!(T, summation) meanAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) centeredSumOfSquares;
Summator!(T, summation) centeredSumOfCubes;
Summator!(T, summation) centeredSumOfQuarts;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
void put()(KurtosisAccumulator!(T, kurtosisAlgo, summation) v);
@property F kurtosis(F = T)(bool isPopulation, bool isRaw)
if (isFloatingPoint!F);
struct KurtosisAccumulator(T, KurtosisAlgo kurtosisAlgo, Summation summation) if (isMutable!T && (kurtosisAlgo == KurtosisAlgo.twoPass || kurtosisAlgo == KurtosisAlgo.threePass));
Examples:
twoPass & threePass
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;
enum RawTrueCT = true;
enum RawFalseCT = false;
bool RawTrueRT = true;
bool RawFalseRT = false;

auto v = KurtosisAccumulator!(double, KurtosisAlgo.twoPass, Summation.naive)(x);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual(38.062853 / 12));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual(38.062853 / 12));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual(38.062853 / 12 - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual(38.062853 / 12 - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));

auto w = KurtosisAccumulator!(double, KurtosisAlgo.threePass, Summation.naive)(x);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual(38.062853 / 12));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual(38.062853 / 12));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual(38.062853 / 12 - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual(38.062853 / 12 - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
this(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);
this(U)(U[] array);
this(T)(T withAsSlice)
if (hasAsSlice!T);
size_t count;
Summator!(T, summation) scaledSumOfQuarts;
@property F kurtosis(F = T)(bool isPopulation, bool isRaw)
if (isFloatingPoint!F);
struct KurtosisAccumulator(T, KurtosisAlgo kurtosisAlgo, Summation summation) if (isMutable!T && (kurtosisAlgo == KurtosisAlgo.assumeZeroMean));
Examples:
assumeZeroMean
import mir.math.common: approxEqual, pow;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;
auto x = a.center;

enum PopulationTrueCT = true;
enum PopulationFalseCT = false;
bool PopulationTrueRT = true;
bool PopulationFalseRT = false;
enum RawTrueCT = true;
enum RawFalseCT = false;
bool RawTrueRT = true;
bool RawFalseRT = false;

KurtosisAccumulator!(double, KurtosisAlgo.assumeZeroMean, Summation.naive) v;
v.put(x);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual(792.784119 / pow(54.765625 / 11, 2.0) * (12.0 * 13.0) / (11.0 * 10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual(792.784119 / pow(54.765625 / 11, 2.0) * (12.0 * 13.0) / (11.0 * 10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual(792.784119 / pow(54.765625 / 11, 2.0) * (12.0 * 13.0) / (11.0 * 10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual(792.784119 / pow(54.765625 / 11, 2.0) * (12.0 * 13.0) / (11.0 * 10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0) + 3.0));

v.put(4.0);
assert(v.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual((1048.784119 / 13) / pow(70.765625 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual((1048.784119 / 13) / pow(70.765625 / 13, 2.0)));
assert(v.kurtosis(PopulationTrueRT, RawFalseRT).approxEqual((1048.784119 / 13) / pow(70.765625 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationTrueCT, RawFalseCT).approxEqual((1048.784119 / 13) / pow(70.765625 / 13, 2.0) - 3.0));
assert(v.kurtosis(PopulationFalseRT, RawFalseRT).approxEqual(1048.784119 / pow(70.765625 / 12, 2.0) * (13.0 * 14.0) / (12.0 * 11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseCT, RawFalseCT).approxEqual(1048.784119 / pow(70.765625 / 12, 2.0) * (13.0 * 14.0) / (12.0 * 11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0)));
assert(v.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual(1048.784119 / pow(70.765625 / 12, 2.0) * (13.0 * 14.0) / (12.0 * 11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
assert(v.kurtosis(PopulationFalseCT, RawTrueCT).approxEqual(1048.784119 / pow(70.765625 / 12, 2.0) * (13.0 * 14.0) / (12.0 * 11.0 * 10.0) - 3.0 * (12.0 * 12.0) / (11.0 * 10.0) + 3.0));
this(Range)(Range r)
if (isIterable!Range);
this()(T x);
VarianceAccumulator!(T, VarianceAlgo.assumeZeroMean, summation) varianceAccumulator;
@property size_t count();
@property F mean(F = T)();
Summator!(T, summation) centeredSumOfQuarts;
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
void put()(KurtosisAccumulator!(T, kurtosisAlgo, summation) v);
@property F kurtosis(F = T)(bool isPopulation, bool isRaw)
if (isFloatingPoint!F);
template kurtosis(F, KurtosisAlgo kurtosisAlgo = KurtosisAlgo.online, Summation summation = Summation.appropriate)

template kurtosis(KurtosisAlgo kurtosisAlgo = KurtosisAlgo.online, Summation summation = Summation.appropriate)

template kurtosis(F, string kurtosisAlgo, string summation = "appropriate")

template kurtosis(string kurtosisAlgo, string summation = "appropriate")
Calculates the kurtosis of the input
By default, if F is not floating point type, then the result will have a double type if F is implicitly convertible to a floating point type.
Parameters:
F controls type of output
kurtosisAlgo algorithm for calculating kurtosis (default: KurtosisAlgo.online) summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The kurtosis of the input, must be floating point or complex type
Examples:
Simple example
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

assert(kurtosis([1.0, 2, 3, 4]).approxEqual(-1.2));

assert(kurtosis([1.0, 2, 4, 5]).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) * (3.0 * 5.0) / (2.0 * 1.0) - 3.0 * (3.0 * 3.0) / (2.0 * 1.0)));
assert(kurtosis([1.0, 2, 4, 5], true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) - 3.0));
assert(kurtosis([1.0, 2, 4, 5], false, true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0) * (3.0 * 5.0) / (2.0 * 1.0) - 3.0 * (3.0 * 3.0) / (2.0 * 1.0) + 3.0));
assert(kurtosis([1.0, 2, 4, 5], true, true).approxEqual((34.0 / 4) / pow(10.0 / 4, 2.0)));

assert(kurtosis!float([0, 1, 2, 3, 4, 6].sliced(3, 2)).approxEqual(-0.2999999));

static assert(is(typeof(kurtosis!float([1, 2, 3])) == float));
Examples:
Kurtosis of vector
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.kurtosis.approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
Examples:
Kurtosis of matrix
import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.kurtosis.approxEqual((792.784119 / 12) / pow(54.765625 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));
Examples:
Column kurtosis of matrix
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual, pow;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: alongDim, byDim, map;

auto x = [
    [0.0,  1.0,  1.5, 2.0], 
    [3.5, 4.25,  2.0, 7.5],
    [5.0,  1.0,  1.5, 0.0],
    [1.5,  4.5, 4.75, 0.5]
].fuse;
auto result = [-2.067182, -5.918089, 3.504056, 2.690240];

// Use byDim or alongDim with map to compute kurtosis of row/column.
assert(x.byDim!1.map!kurtosis.all!approxEqual(result));
assert(x.alongDim!0.map!kurtosis.all!approxEqual(result));

// FIXME
// Without using map, computes the kurtosis of the whole slice
// assert(x.byDim!1.kurtosis == x.sliced.kurtosis);
// assert(x.alongDim!0.kurtosis == x.sliced.kurtosis);
Examples:
Can also set algorithm type
import mir.math.common: approxEqual, pow;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto x = a + 100_000_000_000;

// The default online algorithm is numerically unstable in this case
auto y = x.kurtosis;
assert(!y.approxEqual((792.78411865 / 12) / pow(54.76562500 / 12, 2.0) * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)));

// The naive algorithm has an assert error in this case because standard
// deviation is calculated naively as zero. The kurtosis formula would then
// be dividing by zero. 
//auto z0 = x.kurtosis!(real, "naive");

// The two-pass algorithm is also numerically unstable in this case
auto z1 = x.kurtosis!"twoPass";
assert(!z1.approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(!z1.approxEqual(y));

// However, the three-pass algorithm is numerically stable in this case
auto z2 = x.kurtosis!"threePass";
assert(z2.approxEqual(38.062853 / 12 * (11.0 * 13.0) / (10.0 * 9.0) - 3.0 * (11.0 * 11.0) / (10.0 * 9.0)) + 3.0);
assert(!z2.approxEqual(y));

// And the assumeZeroMean algorithm provides the incorrect answer, as expected
auto z3 = x.kurtosis!"assumeZeroMean";
assert(!z3.approxEqual(y));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

// Set population/sample kurtosis, excess/raw kurtosis, kurtosis algorithm,
// sum algorithm or output type

auto a = [1.0, 1e72, 1, -1e72].sliced;
auto x = a * 10_000;

bool PopulationTrueRT = true;
bool PopulationFalseRT = false;
enum PopulationTrueCT = true;

enum RawTrueCT = true;
bool RawTrueRT = true;
bool RawFalseRT = false;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and taking `x` to the fourth power, the first and third numbers in the slice
have precision too low to be included in the centered sum of cubes. 
+/
assert(x.kurtosis.approxEqual(1.5));
assert(x.kurtosis(PopulationFalseRT).approxEqual(1.5));
assert(x.kurtosis(PopulationTrueRT).approxEqual(-1.0));
assert(x.kurtosis(PopulationTrueCT).approxEqual(-1.0));
assert(x.kurtosis(PopulationTrueRT, RawTrueRT).approxEqual(2.0));
assert(x.kurtosis(PopulationFalseRT, RawTrueRT).approxEqual(4.5));
assert(x.kurtosis(PopulationTrueCT, RawTrueCT).approxEqual(2.0));

assert(x.kurtosis!("online").approxEqual(1.5));
assert(x.kurtosis!("online", "kbn").approxEqual(1.5));
assert(x.kurtosis!("online", "kb2").approxEqual(1.5));
assert(x.kurtosis!("online", "precise").approxEqual(1.5));
assert(x.kurtosis!(double, "online", "precise").approxEqual(1.5));
assert(x.kurtosis!(double, "online", "precise")(PopulationTrueRT).approxEqual(-1.0));
assert(x.kurtosis!(double, "online", "precise")(PopulationTrueRT, RawTrueRT).approxEqual(2.0));

auto y = [uint.max - 3, uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.kurtosis!(ulong, "threePass");
assert(z.approxEqual(-1.2));
static assert(is(typeof(z) == double));
Examples:
For integral slices, can pass output type as template parameter to ensure output type is correct.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.kurtosis;
assert(y.approxEqual(0.223394));
static assert(is(typeof(y) == double));

assert(x.kurtosis!float.approxEqual(0.223394));
Examples:
Kurtosis works for other user-defined types (provided they can be converted to a floating point)
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

static struct Foo {
    float x;
    alias x this;
}

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f), Foo(4f)];
assert(foo.kurtosis.approxEqual(-1.2f));
Examples:
Compute kurtosis along specified dimention of tensors
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: as, iota, alongDim, map, repeat;

auto x = [
    [0.0,  1,  3,  5],
    [3.0,  4,  5,  7],
    [6.0,  7, 10, 11],
    [9.0, 12, 15, 12]
].fuse;

assert(x.kurtosis.approxEqual(-0.770040));

auto m0 = [-1.200000, -0.152893, -1.713859, -3.869005];
assert(x.alongDim!0.map!kurtosis.all!approxEqual(m0));
assert(x.alongDim!(-2).map!kurtosis.all!approxEqual(m0));

auto m1 = [-1.699512, 0.342857, -4.339100, 1.500000];
assert(x.alongDim!1.map!kurtosis.all!approxEqual(m1));
assert(x.alongDim!(-1).map!kurtosis.all!approxEqual(m1));

assert(iota(4, 5, 6, 7).as!double.alongDim!0.map!kurtosis.all!approxEqual(repeat(-1.2, 5, 6, 7)));
Examples:
Arbitrary kurtosis
import mir.math.common: approxEqual;

assert(kurtosis(1.0, 2, 3, 4).approxEqual(-1.2));
assert(kurtosis!float(1, 2, 3, 4).approxEqual(-1.2f));
stdevType!F kurtosis(Range)(Range r, bool isPopulation = false, bool isRaw = false)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
bool isPopulation true if population kurtosis, false if sample kurtosis (default)
bool isRaw true if raw kurtosis, false if excess kurtosis (default)
stdevType!F kurtosis(scope const F[] ar...);
Parameters:
F[] ar values
struct EntropyAccumulator(T, Summation summation);
Summator!(T, summation) summator;
const pure nothrow @nogc @property @safe F entropy(F = T)();
void put(Range)(Range r)
if (isIterable!Range);
void put()(T x);
void put(U)(EntropyAccumulator!(U, summation) e);
template entropyType(T)
template entropy(F, Summation summation = Summation.appropriate)
Computes the entropy of the input.
By default, if F is not a floating point type, then the result will have a double type if F is implicitly convertible to a floating point type.
Parameters:
F controls type of output
summation algorithm for summing the individual entropy values (default: Summation.appropriate)
Returns:
The entropy of all the elements in the input, must be floating point type
See Also:
entropyType!Range entropy(Range)(Range r)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
entropyType!F entropy(scope const F[] ar...);
Parameters:
F[] ar values
template entropy(Summation summation = Summation.appropriate)

template entropy(F, string summation)

template entropy(string summation)
Examples:
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(entropy([0.166667, 0.333333, 0.50]).approxEqual(-1.011404));

assert(entropy!float([0.05, 0.1, 0.15, 0.2, 0.25, 0.25].sliced(3, 2)).approxEqual(-1.679648));

static assert(is(typeof(entropy!float([0.166667, 0.333333, 0.50])) == float));
Examples:
Entropy of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

double[] a = [1.0, 2, 3,  4,  5,  6, 7, 8, 9, 10, 11, 12];
a[] /= 78.0;

auto x = a.sliced;
assert(x.entropy.approxEqual(-2.327497));
Examples:
Entropy of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

double[] a = [1.0, 2, 3,  4,  5,  6, 7, 8, 9, 10, 11, 12];
a[] /= 78.0;

auto x = a.fuse;
assert(x.entropy.approxEqual(-2.327497));
Examples:
Column entropy of matrix
import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.topology: alongDim, byDim, map;

double[][] a = [
    [1.0, 2, 3,  4,  5,  6], 
    [7.0, 8, 9, 10, 11, 12]
];
a[0][] /= 78.0;
a[1][] /= 78.0;

auto x = a.fuse;
auto result = [-0.272209, -0.327503, -0.374483, -0.415678, -0.452350, -0.485273];

// Use byDim or alongDim with map to compute entropy of row/column.
assert(x.byDim!1.map!entropy.all!approxEqual(result));
assert(x.alongDim!0.map!entropy.all!approxEqual(result));

// FIXME
// Without using map, computes the entropy of the whole slice
// assert(x.byDim!1.entropy == x.sliced.entropy);
// assert(x.alongDim!0.entropy == x.sliced.entropy);
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

auto a = [1, 1e100, 1, 1e100].sliced;

auto x = a * 10_000;

assert(x.entropy!"kbn".approxEqual(4.789377e106));
assert(x.entropy!"kb2".approxEqual(4.789377e106));
assert(x.entropy!"precise".approxEqual(4.789377e106));
assert(x.entropy!(double, "precise").approxEqual(4.789377e106));
Examples:
For integral slices, pass output type as template parameter to ensure output type is correct.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [3, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 3].sliced;

auto y = x.entropy;
assert(y.approxEqual(43.509472));
static assert(is(typeof(y) == double));

assert(x.entropy!float.approxEqual(43.509472f));
Examples:
Arbitrary entropy
import mir.math.common: approxEqual;

assert(entropy(0.25, 0.25, 0.25, 0.25).approxEqual(-1.386294));
assert(entropy!float(0.25, 0.25, 0.25, 0.25).approxEqual(-1.386294));
entropyType!Range entropy(Range)(Range r)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
entropyType!T entropy(T)(scope const T[] ar...);
Parameters:
T[] ar values
template coefficientOfVariation(F, VarianceAlgo varianceAlgo = VarianceAlgo.online, Summation summation = Summation.appropriate)

template coefficientOfVariation(VarianceAlgo varianceAlgo = VarianceAlgo.online, Summation summation = Summation.appropriate)
Calculates the coefficient of variation of the input.
The coefficient of variation is calculated by dividing either the population or sample (default) standard deviation by the mean of the input. According to wikipedia, "the coefficient of variation should be computed computed for data measured on a ratio scale, that is, scales that have a meaningful zero and hence allow for relative comparison of two measurements." In addition, for "small- and moderately-sized datasets", the coefficient of variation is biased, even when using the sample standard deviation.
By default, if F is not floating point type, then the result will have a double type if F is implicitly convertible to a floating point type.
Parameters:
F controls type of output
varianceAlgo algorithm for calculating variance (default: VarianceAlgo.online) summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The coefficient of varition of the input, must be floating point type
See Also: Coefficient of variation.
stdevType!F coefficientOfVariation(Range)(Range r, bool isPopulation = false)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
bool isPopulation true if population variance, false if sample variance (default)
stdevType!F coefficientOfVariation(scope const F[] ar...);
Parameters:
F[] ar values
template coefficientOfVariation(F, string varianceAlgo, string summation = "appropriate")

template coefficientOfVariation(string varianceAlgo, string summation = "appropriate")
Examples:
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(coefficientOfVariation([1.0, 2, 3]).approxEqual(1.0 / 2.0));
assert(coefficientOfVariation([1.0, 2, 3], true).approxEqual(0.816497 / 2.0));

assert(coefficientOfVariation!float([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(1.870829 / 2.5));

static assert(is(typeof(coefficientOfVariation!float([1, 2, 3])) == float));
Examples:
Coefficient of variation of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.coefficientOfVariation.approxEqual(2.231299 / 2.437500));
Examples:
Coefficient of variation of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.coefficientOfVariation.approxEqual(2.231299 / 2.437500));
Examples:
Can also set algorithm type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto x = a + 1_000_000_000;

auto y = x.coefficientOfVariation;
assert(y.approxEqual(2.231299 / 1_000_000_002.437500));

// The naive variance algorithm is numerically unstable in this case, but
// the difference is small as coefficientOfVariation is a ratio
auto z0 = x.coefficientOfVariation!"naive";
assert(!z0.approxEqual(y, 0x1p-20f, 0x1p-30f));

// But the two-pass algorithm provides a consistent answer
auto z1 = x.coefficientOfVariation!"twoPass";
assert(z1.approxEqual(y));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

// Set population standard deviation, standardDeviation algorithm, sum algorithm or output type

auto a = [1.0, 1e100, 1, -1e100].sliced;
auto x = a * 10_000;

bool populationTrue = true;

/++
For this case, failing to use a summation algorithm results in an assert
error because the mean is zero due to floating point precision issues.
+/
//assert(x.coefficientOfVariation!("online").approxEqual(8.164966e103 / 0.0));

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.coefficientOfVariation!("online", "kbn").approxEqual(8.164966e103 / 5000.0));
assert(x.coefficientOfVariation!("online", "kb2").approxEqual(8.164966e103 / 5000.0));
assert(x.coefficientOfVariation!("online", "precise").approxEqual(8.164966e103 / 5000.0));
assert(x.coefficientOfVariation!(double, "online", "precise").approxEqual(8.164966e103 / 5000.0));
assert(x.coefficientOfVariation!(double, "online", "precise")(populationTrue).approxEqual(7.071068e103 / 5000.0));


auto y = [uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.coefficientOfVariation!ulong;
assert(z == (1.0 / (cast(double) uint.max - 1)));
static assert(is(typeof(z) == double));
assert(y.coefficientOfVariation!(ulong, "online") == (1.0 / (cast(double) uint.max - 1)));
Examples:
For integral slices, pass output type as template parameter to ensure output type is correct.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.coefficientOfVariation;
assert(y.approxEqual(2.151462f / 2.416667));
static assert(is(typeof(y) == double));

assert(x.coefficientOfVariation!float.approxEqual(2.151462f / 2.416667));
Examples:
coefficientOfVariation works for other user-defined types (provided they can be converted to a floating point)
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

static struct Foo {
    float x;
    alias x this;
}

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f)];
assert(foo.coefficientOfVariation.approxEqual(1f / 2f));
Examples:
Arbitrary coefficientOfVariation
import mir.math.common: approxEqual;

assert(coefficientOfVariation(1.0, 2, 3).approxEqual(1.0 / 2.0));
assert(coefficientOfVariation!float(1, 2, 3).approxEqual(1f / 2f));
struct MomentAccumulator(T, size_t N, Summation summation) if (N > 0 && isMutable!T);
Examples:
Raw moment
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;
auto x = a.center;

MomentAccumulator!(double, 2, Summation.naive) v;
v.put(x);

assert(v.moment.approxEqual(54.76562 / 12));

v.put(4.0);
assert(v.moment.approxEqual(70.76562 / 13));
Examples:
Central moment
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

MomentAccumulator!(double, 2, Summation.naive) v;
auto m = mean(x);
v.put(x, m);
assert(v.moment.approxEqual(54.76562 / 12));
Examples:
Standardized moment with scaled calculation
import mir.math.common: approxEqual, sqrt;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto u = VarianceAccumulator!(double, VarianceAlgo.twoPass, Summation.naive)(x);
MomentAccumulator!(double, 3, Summation.naive) v;
v.put(x, u.mean, u.variance(true).sqrt);
//assert(v.moment.approxEqual(12.000999 / 12));
assert(v.count == 12);
Summator!(T, summation) summator;
size_t count;
const pure nothrow @nogc @property @safe F moment(F = T)();
const pure nothrow @nogc @property @safe F sumOfPower(F = T)();
void put(Range)(Range r)
if (isIterable!Range);
void put(Range)(Range r, T m)
if (isIterable!Range);
void put(Range)(Range r, T m, T s)
if (isIterable!Range);
void put()(T x);
void put()(MomentAccumulator!(T, N, summation) m);
this(Range)(Range r)
if (isIterable!Range);
this(Range)(Range r, T m)
if (isIterable!Range);
this(Range)(Range r, T m, T s)
if (isIterable!Range);
this()(T x);
this()(T x, T m);
this()(T x, T m, T s);
template rawMoment(F, size_t N, Summation summation = Summation.appropriate) if (N > 0)

template rawMoment(size_t N, Summation summation = Summation.appropriate) if (N > 0)

template rawMoment(F, size_t N, string summation) if (N > 0)

template rawMoment(size_t N, string summation) if (N > 0)
Calculates the n-th raw moment of the input.
By default, if F is not floating point type or complex type, then the result will have a double type if F is implicitly convertible to a floating point type or have a cdouble type if F is implicitly convertible to a complex type.
Parameters:
F controls type of output
N controls n-th raw moment summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The n-th raw moment of the input, must be floating point or complex type
Examples:
Basic implementation
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(rawMoment!2([1.0, 2, 3]).approxEqual(14.0 / 3));
assert(rawMoment!3([1.0, 2, 3]).approxEqual(36.0 / 3));

assert(rawMoment!(float, 2)([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(55f / 6));
static assert(is(typeof(rawMoment!(float, 2)([1, 2, 3])) == float));
Examples:
Raw Moment of vector
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;
auto x = a.center;

assert(x.rawMoment!2.approxEqual(54.76562 / 12));
Examples:
Raw Moment of matrix
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.fuse: fuse;

auto a = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;
auto x = a.center;

assert(x.rawMoment!2.approxEqual(54.76562 / 12));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

//Set sum algorithm or output type

auto a = [1.0, 1e100, 1, -1e100].sliced;
auto b = a * 10_000;
auto x = b.center;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.rawMoment!2.approxEqual(2.0e208 / 4));

assert(x.rawMoment!(2, "kbn").approxEqual(2.0e208 / 4));
assert(x.rawMoment!(2, "kb2").approxEqual(2.0e208 / 4));
assert(x.rawMoment!(2, "precise").approxEqual(2.0e208 / 4));
assert(x.rawMoment!(double, 2, "precise").approxEqual(2.0e208 / 4));

auto y = uint.max.repeat(3);
auto z = y.rawMoment!(ulong, 2);
assert(z.approxEqual(cast(double) (cast(ulong) uint.max) ^^ 2u));
static assert(is(typeof(z) == double));
Examples:
rawMoment works for complex numbers and other user-defined types (provided they can be converted to a floating point or complex type)
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [1.0 + 2i, 2 + 3i, 3 + 4i, 4 + 5i].sliced;
assert(x.rawMoment!2.approxEqual((-24 + 80.0i)/ 4));
Examples:
Arbitrary raw moment
import mir.math.common: approxEqual;

assert(rawMoment!2(1.0, 2, 3).approxEqual(14.0 / 3));
assert(rawMoment!(float, 2)(1, 2, 3).approxEqual(14f / 3));
meanType!F rawMoment(Range)(Range r)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
meanType!F rawMoment(scope const F[] ar...);
Parameters:
F[] ar values
template centralMoment(F, size_t N, Summation summation = Summation.appropriate) if (N > 0)

template centralMoment(size_t N, Summation summation = Summation.appropriate) if (N > 0)

template centralMoment(F, size_t N, string summation) if (N > 0)

template centralMoment(size_t N, string summation) if (N > 0)
Calculates the n-th central moment of the input.
By default, if F is not floating point type or complex type, then the result will have a double type if F is implicitly convertible to a floating point type or have a cdouble type if F is implicitly convertible to a complex type.
Parameters:
F controls type of output
N controls n-th central moment summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The n-th central moment of the input, must be floating point or complex type
Examples:
Basic implementation
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(centralMoment!2([1.0, 2, 3]).approxEqual(2.0 / 3));
assert(centralMoment!3([1.0, 2, 3]).approxEqual(0.0 / 3));

assert(centralMoment!(float, 2)([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(17.5f / 6));
static assert(is(typeof(centralMoment!(float, 2)([1, 2, 3])) == float));
Examples:
Central Moment of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.centralMoment!2.approxEqual(54.76562 / 12));
Examples:
Central Moment of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.centralMoment!2.approxEqual(54.76562 / 12));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.math.stat: center;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

//Set sum algorithm or output type

auto a = [1.0, 1e100, 1, -1e100].sliced;
auto b = a * 10_000;
auto x = b.center;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.centralMoment!2.approxEqual(2.0e208 / 4));

assert(x.centralMoment!(2, "kbn").approxEqual(2.0e208 / 4));
assert(x.centralMoment!(2, "kb2").approxEqual(2.0e208 / 4));
assert(x.centralMoment!(2, "precise").approxEqual(2.0e208 / 4));
assert(x.centralMoment!(double, 2, "precise").approxEqual(2.0e208 / 4));

auto y = uint.max.repeat(3);
auto z = y.centralMoment!(ulong, 2);
assert(z.approxEqual(0.0));
static assert(is(typeof(z) == double));
Examples:
centralMoment works for complex numbers and other user-defined types (provided they can be converted to a floating point or complex type)
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [1.0 + 2i, 2 + 3i, 3 + 4i, 4 + 5i].sliced;
assert(x.centralMoment!2.approxEqual((0.0 + 10.0i) / 4));
Examples:
Arbitrary central moment
import mir.math.common: approxEqual;

assert(centralMoment!2(1.0, 2, 3).approxEqual(2.0 / 3));
assert(centralMoment!(float, 2)(1, 2, 3).approxEqual(2f / 3));
meanType!F centralMoment(Range)(Range r)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
meanType!F centralMoment(scope const F[] ar...);
Parameters:
F[] ar values
enum StandardizedMomentAlgo: int;
scaled
Calculates n-th standardized moment as E(((x - u) / sigma) ^^ N)
centered
Calculates n-th standardized moment as E(((x - u) ^^ N) / ((x - u) ^^ (N / 2)))
template standardizedMoment(F, size_t N, StandardizedMomentAlgo standardizedMomentAlgo = StandardizedMomentAlgo.scaled, VarianceAlgo varianceAlgo = VarianceAlgo.twoPass, Summation summation = Summation.appropriate) if (N > 0)

template standardizedMoment(size_t N, StandardizedMomentAlgo standardizedMomentAlgo = StandardizedMomentAlgo.scaled, VarianceAlgo varianceAlgo = VarianceAlgo.twoPass, Summation summation = Summation.appropriate) if (N > 0)

template standardizedMoment(F, size_t N, string standardizedMomentAlgo, string varianceAlgo = "twoPass", string summation = "appropriate") if (N > 0)

template standardizedMoment(size_t N, string standardizedMomentAlgo, string varianceAlgo = "twoPass", string summation = "appropriate") if (N > 0)
Calculates the n-th standardized moment of the input.
By default, if F is not floating point type or complex type, then the result will have a double type if F is implicitly convertible to a floating point type or have a cdouble type if F is implicitly convertible to a complex type.
Parameters:
F controls type of output
N controls n-th standardized moment summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The n-th standardized moment of the input, must be floating point or complex type
Examples:
Basic implementation
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(standardizedMoment!1([1.0, 2, 3]).approxEqual(0.0));
assert(standardizedMoment!2([1.0, 2, 3]).approxEqual(1.0));
assert(standardizedMoment!3([1.0, 2, 3]).approxEqual(0.0 / 3));
assert(standardizedMoment!4([1.0, 2, 3]).approxEqual(4.5 / 3));

assert(standardizedMoment!(float, 2)([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(6f / 6));
static assert(is(typeof(standardizedMoment!(float, 2)([1, 2, 3])) == float));
Examples:
Standardized Moment of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.standardizedMoment!3.approxEqual(12.000999 / 12));
Examples:
Standardized Moment of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.standardizedMoment!3.approxEqual(12.000999 / 12));
Examples:
Can also set algorithm type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto a = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

auto x = a + 100_000_000_000;

// The default algorithm is numerically stable in this case
auto y = x.standardizedMoment!3;
assert(y.approxEqual(12.000999 / 12));

// The online algorithm is numerically unstable in this case
auto z1 = x.standardizedMoment!(3, "scaled", "online");
assert(!z1.approxEqual(12.000999 / 12));
assert(!z1.approxEqual(y));

// It is also numerically unstable when using StandardizedMomentAlgo.centered
auto z2 = x.standardizedMoment!(3, "centered", "online");
assert(!z2.approxEqual(12.000999 / 12));
assert(!z2.approxEqual(y));
Examples:
Can also set algorithm or output type
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

//Set standardized moment algorithm, variance algorithm, sum algorithm, or output type

auto a = [1.0, 1e98, 1, -1e98].sliced;
auto x = a * 10_000;

/++
Due to Floating Point precision, when centering `x`, subtracting the mean 
from the second and fourth numbers has no effect. Further, after centering 
and squaring `x`, the first and third numbers in the slice have precision 
too low to be included in the centered sum of squares. 
+/
assert(x.standardizedMoment!3.approxEqual(0.0));

assert(x.standardizedMoment!(3, "scaled", "online").approxEqual(0.0));
assert(x.standardizedMoment!(3, "centered", "online").approxEqual(0.0));
assert(x.standardizedMoment!(3, "scaled", "online", "kbn").approxEqual(0.0));
assert(x.standardizedMoment!(3, "scaled", "online", "kb2").approxEqual(0.0));
assert(x.standardizedMoment!(3, "scaled", "online", "precise").approxEqual(0.0));
assert(x.standardizedMoment!(double, 3, "scaled", "online", "precise").approxEqual(0.0));

auto y = [uint.max - 2, uint.max - 1, uint.max].sliced;
auto z = y.standardizedMoment!(ulong, 3);
assert(z == 0.0);
static assert(is(typeof(z) == double));
Examples:
For integral slices, can pass output type as template parameter to ensure output type is correct. By default, they get converted to double.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.standardizedMoment!3;
assert(y.approxEqual(9.666455 / 12));
static assert(is(typeof(y) == double));

assert(x.standardizedMoment!(float, 3).approxEqual(9.666455f / 12));
Examples:
Arbitrary standardized moment
import mir.math.common: approxEqual;

assert(standardizedMoment!3(1.0, 2, 3).approxEqual(0.0 / 3));
assert(standardizedMoment!(float, 3)(1, 2, 3).approxEqual(0f / 3));
assert(standardizedMoment!(float, 3, "centered")(1, 2, 3).approxEqual(0f / 3));
stdevType!F standardizedMoment(Range)(Range r)
if (isIterable!Range);
Parameters:
Range r range, must be finite iterable
stdevType!F standardizedMoment(scope const F[] ar...);
Parameters:
F[] ar values
enum MomentAlgo: int;
raw
nth raw moment, E(x ^^ n)
central
nth central moment, E((x - u) ^^ n)
standardized
nth standardized moment, E(((x - u) / sigma) ^^ n)
template moment(F, size_t N, MomentAlgo momentAlgo, Summation summation = Summation.appropriate)

template moment(size_t N, MomentAlgo momentAlgo, Summation summation = Summation.appropriate)

template moment(F, size_t N, string momentAlgo, string summation = "appropriate")

template moment(size_t N, string momentAlgo, string summation = "appropriate")
Calculates the n-th moment of the input.
Parameters:
F controls type of output
N controls n-th standardized moment momentAlgo: type of moment to be calculated summation: algorithm for calculating sums (default: Summation.appropriate)
Returns:
The n-th moment of the input, must be floating point or complex type
Examples:
Basic implementation
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

assert(moment!(1, "raw")([1.0, 2, 3]).approxEqual(6.0 / 3));
assert(moment!(2, "raw")([1.0, 2, 3]).approxEqual(14.0 / 3));
assert(moment!(3, "raw")([1.0, 2, 3]).approxEqual(36.0 / 3));
assert(moment!(4, "raw")([1.0, 2, 3]).approxEqual(98.0 / 3));

assert(moment!(1, "central")([1.0, 2, 3]).approxEqual(0.0 / 3));
assert(moment!(2, "central")([1.0, 2, 3]).approxEqual(2.0 / 3));
assert(moment!(3, "central")([1.0, 2, 3]).approxEqual(0.0 / 3));
assert(moment!(4, "central")([1.0, 2, 3]).approxEqual(2.0 / 3));

assert(moment!(1, "standardized")([1.0, 2, 3]).approxEqual(0.0));
assert(moment!(2, "standardized")([1.0, 2, 3]).approxEqual(1.0));
assert(moment!(3, "standardized")([1.0, 2, 3]).approxEqual(0.0 / 3));
assert(moment!(4, "standardized")([1.0, 2, 3]).approxEqual(4.5 / 3));

assert(moment!(float, 2, "standardized")([0, 1, 2, 3, 4, 5].sliced(3, 2)).approxEqual(6f / 6));
static assert(is(typeof(moment!(float, 2, "standardized")([1, 2, 3])) == float));
Examples:
Standardized Moment of vector
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0.0, 1.0, 1.5, 2.0, 3.5, 4.25,
          2.0, 7.5, 5.0, 1.0, 1.5, 0.0].sliced;

assert(x.moment!(3, "standardized").approxEqual(12.000999 / 12));
Examples:
Standardized Moment of matrix
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;

auto x = [
    [0.0, 1.0, 1.5, 2.0, 3.5, 4.25],
    [2.0, 7.5, 5.0, 1.0, 1.5, 0.0]
].fuse;

assert(x.moment!(3, "standardized").approxEqual(12.000999 / 12));
Examples:
For integral slices, can pass output type as template parameter to ensure output type is correct. By default, they get converted to double.
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

auto x = [0, 1, 1, 2, 4, 4,
          2, 7, 5, 1, 2, 0].sliced;

auto y = x.moment!(3, "standardized");
assert(y.approxEqual(9.666455 / 12));
static assert(is(typeof(y) == double));

assert(x.moment!(float, 3, "standardized").approxEqual(9.666455f / 12));
Examples:
Arbitrary standardized moment
import mir.math.common: approxEqual;

assert(moment!(3, "standardized")(1.0, 2, 3).approxEqual(0.0 / 3));
assert(moment!(float, 3, "standardized")(1, 2, 3).approxEqual(0f / 3));
meanType!F moment(Range)(Range r)
if (isIterable!Range && (momentAlgo != MomentAlgo.standardized));
Parameters:
Range r range, must be finite iterable
stdevType!F moment(Range)(Range r)
if (isIterable!Range && (momentAlgo == MomentAlgo.standardized));
Parameters:
Range r range, must be finite iterable
meanType!F moment()(scope const F[] ar...)
if (momentAlgo != MomentAlgo.standardized);
Parameters:
F[] ar values
stdevType!F moment()(scope const F[] ar...)
if (momentAlgo == MomentAlgo.standardized);
Parameters:
F[] ar values