If you spot a problem with this page, click here to create a GitHub issue.

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.univariate

This module contains algorithms for univariate descriptive statistics.


Category	Symbols
Location	*gmean* *hmean* *mean* *median*
Deviation	dispersion entropy interquartileRange medianAbsoluteDeviation quantile *standardDeviation* *variance*
Higher Moments, etc.	kurtosis skewness
Other Moment Functions	centralMoment coefficientOfVariation moment rawMoment standardizedMoment
Accumulators	EntropyAccumulator *GMeanAccumulator* KurtosisAccumulator *MeanAccumulator* MomentAccumulator SkewnessAccumulator *VarianceAccumulator*
Algorithms	KurtosisAlgo MomentAlgo QuantileAlgo SkewnessAlgo StandardizedMomentAlgo *VarianceAlgo*
Types	entropyType *gmeanType* *meanType* quantileType *statType* *stdevType*

License:

Apache-2.0

Authors:

John Michael Hall, Ilya Yaroshenko

public import mir.math.stat : gmean, GMeanAccumulator, hmean, mean, MeanAccumulator, meanType, median, standardDeviation, statType, stdevType, variance, VarianceAccumulator, VarianceAlgo;

public import mir.math.sum : Summation;

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.


Type	m	gamma
*Discontinuous sample quantile*
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*
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:

quantile

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 quantileType(T, QuantileAlgo quantileAlgo)

For all QuantileAlgo except QuantileAlgo.type1 and QuantileAlgo.type3, this is an alias to the meanType of T

For QuantileAlgo.type1 and QuantileAlgo.type3, this is an alias to the elementType of T.

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 compile-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 a type for which isComplex!F is true.

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:

median , partitionAt , elementType

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)

Examples:

Simple example

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;
          
assert(x.quantile(0.5).approxEqual(2.0));

auto qtile = [0.25, 0.75].sliced;

assert(x.quantile(qtile).all!approxEqual([1.0, 3.0]));
assert(x.quantile(0.25, 0.75).all!approxEqual([1.0, 3.0]));

Examples:

Modify probability in place

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 qtile = [0.25, 0.75].sliced;
auto qtile_copy = qtile.dup;

x.quantile!("type7", false, true)(qtile);
assert(qtile.all!approxEqual([1.0, 3.0]));
assert(!qtile.all!approxEqual(qtile_copy));

Examples:

Quantile of vector

import mir.algorithm.iteration: all;
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.quantile(0.5).approxEqual(5.20));

auto qtile = [0.25, 0.75].sliced;

assert(x.quantile(qtile).all!approxEqual([3.250, 8.500]));

Examples:

Quantile of matrix

import mir.algorithm.iteration: all;
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.quantile(0.5).approxEqual(5.20));

auto qtile = [0.25, 0.75].sliced;

assert(x.quantile(qtile).all!approxEqual([3.250, 8.500]));

Examples:

Row quantile of matrix

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;
import mir.ndslice.fuse: fuse;
import mir.ndslice.slice: sliced;
import mir.ndslice.topology: alongDim, byDim, map, flattened;

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;

auto result0 = [5.200, 5.700];

// Use byDim or alongDim with map to compute median of row/column.
assert(x.byDim!0.map!(a => a.quantile(0.5)).all!approxEqual(result0));
assert(x.alongDim!1.map!(a => a.quantile(0.5)).all!approxEqual(result0));

auto qtile = [0.25, 0.75].sliced;
auto result1 = [[3.750, 7.600], [2.825, 9.025]];

assert(x.byDim!0.map!(a => a.quantile(qtile)).all!(all!approxEqual)(result1));

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.quantile!(QuantileAlgo.type7, true)(0.5);
assert(!x.all!approxEqual(x_copy));

Examples:

Double-check probability is not modified

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 qtile = [0.25, 0.75].sliced;
auto qtile_copy = qtile.dup;

auto result = x.quantile!("type7", false, false)(qtile);
assert(result.all!approxEqual([1.0, 3.0]));
assert(qtile.all!approxEqual(qtile_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.quantile!"type1"(0.5).approxEqual(5.20));
assert(x.quantile!"type2"(0.5).approxEqual(5.20));
assert(x.quantile!"type3"(0.5).approxEqual(5.20));
assert(x.quantile!"type4"(0.5).approxEqual(5.20));
assert(x.quantile!"type5"(0.5).approxEqual(5.20));
assert(x.quantile!"type6"(0.5).approxEqual(5.20));
assert(x.quantile!"type7"(0.5).approxEqual(5.20));
assert(x.quantile!"type8"(0.5).approxEqual(5.20));
assert(x.quantile!"type9"(0.5).approxEqual(5.20));

Examples:

Can also set algorithm or output type

import mir.ndslice.slice: sliced;

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

auto x = a * 10_000;

auto result0 = x.quantile!float(0.5);
assert(result0 == 10_000f);
static assert(is(typeof(result0) == float));

auto result1 = x.quantile!(float, "type8")(0.5);
assert(result1 == 10_000f);
static assert(is(typeof(result1) == float));

Examples:

Support for integral and user-defined types for type 1 & 3

import mir.ndslice.slice: sliced;
import mir.ndslice.topology: repeat;

auto x = uint.max.repeat(3);
assert(x.quantile!(uint, "type1")(0.5) == uint.max);
assert(x.quantile!(uint, "type3")(0.5) == uint.max);

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

Foo[] foo = [Foo(1f), Foo(2f), Foo(3f)];
assert(foo.quantile!"type1"(0.5) == 2f);
assert(foo.quantile!"type3"(0.5) == 2f);

Examples:

Compute quantile 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],
    [4.0, 5, 7]
].fuse;

assert(x.quantile(0.5).approxEqual(3.5));

auto m0 = [2.0, 3.0, 5.0];
assert(x.alongDim!0.map!(a => a.quantile(0.5)).all!approxEqual(m0));
assert(x.alongDim!(-2).map!(a => a.quantile(0.5)).all!approxEqual(m0));

auto m1 = [1.0, 5.0];
assert(x.alongDim!1.map!(a => a.quantile(0.5)).all!approxEqual(m1));
assert(x.alongDim!(-1).map!(a => a.quantile(0.5)).all!approxEqual(m1));

assert(iota(2, 3, 4, 5).as!double.alongDim!0.map!(a => a.quantile(0.5)).all!approxEqual(iota([3, 4, 5], 3 * 4 * 5 / 2)));

Examples:

Support for array

import mir.algorithm.iteration: all;
import mir.math.common: approxEqual;

double[] x = [3.0, 1.0, 4.0, 2.0, 0.0];
          
assert(x.quantile(0.5).approxEqual(2.0));

double[] qtile = [0.25, 0.75];

assert(x.quantile(qtile).all!approxEqual([1.0, 3.0]));

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.

By default, this function computes the result using quantile, i.e. result = quantile(x, 0.75) - quantile(x, 0.25). There are also overloads for providing a low value, as in result = quantile(x, 1 - low) - quantile(x, low) and both a low and high value, as in result = quantile(x, high) - quantile(x, low).

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:

quantile

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));
assert(x.interquartileRange(0.25).approxEqual(2.0));
assert(x.interquartileRange(0.25, 0.75).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(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice, F lo = 0.25);

Parameters:

Slice!(Iterator, N, kind) `slice`	slice
F `lo`	low value

quantileType!(F, quantileAlgo) interquartileRange(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice, F lo, F hi);

Parameters:

Slice!(Iterator, N, kind) `slice`	slice
F `lo`	low value
F `hi`	high value

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];

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];

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.complex: Complex;
import mir.complex.math: capproxEqual = approxEqual;
import mir.functional: naryFun;
import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

alias C = Complex!double;

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

assert(dispersion([C(1.0, 3), C(2), C(3)]).capproxEqual(C(-4, -6) / 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];

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];
          
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.ndslice.slice: sliced;
import std.complex: Complex;
import std.math.operations: isClose;

auto x = [Complex!double(1, 2), Complex!double(2, 3), Complex!double(3, 4), Complex!double(4, 5)].sliced;
assert(x.dispersion.isClose(Complex!double(0, 10) / 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

alias SkewnessAlgo = KurtosisAlgo;

Skew algorithms (currently an alias to KurtosisAlgo).

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;
import mir.test;

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];

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;
Summator!(T, summation) sumOfCubes;
void put(Range)(Range r) if (isIterable!Range);
void put()(T x);
@property size_t count();
@property F mean(F = T)();
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];

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;
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 size_t count();
@property F mean(F = T)();
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];

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;
T scaledSumOfCubes;
this(Iterator, size_t N, SliceKind kind)(Slice!(Iterator, N, kind) slice);
this(U)(U[] array);
this(T)(T withAsSlice) if (hasAsSlice!T);
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;
Summator!(T, summation) centeredSumOfCubes;
void put(Range)(Range r) if (isIterable!Range);
void put()(T x);
void put()(SkewnessAccumulator!(T, skewnessAlgo, summation) v);
@property size_t count();
@property F mean(F = T)();
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];

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.

online: Similar to Welford's algorithm for updating variance, but adjusted for skewness and 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).
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 and 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 skewness and 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 skewness and 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];

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));

assert(v.skewness(PopulationTrueRT).approxEqual((117.005859 / 12) / pow(54.765625 / 12, 1.5)));

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));

assert(v.skewness(PopulationTrueRT).approxEqual((100.238166 / 13) / pow(57.019231 / 13, 1.5)));

this(Range)(Range r) if (isIterable!Range);
this()(T x);
SkewnessAccumulator!(T, kurtosisAlgo, summation) skewnessAccumulator;
Summator!(T, summation) sumOfQuarts;
void put(Range)(Range r) if (isIterable!Range);
void put()(T x);
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];

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);
SkewnessAccumulator!(T, kurtosisAlgo, summation) skewnessAccumulator;
Summator!(T, summation) centeredSumOfQuarts;
void put(Range)(Range r) if (isIterable!Range);
void put()(T x);
void put()(KurtosisAccumulator!(T, kurtosisAlgo, summation) v);
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];

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.twoPass, Summation.naive) v = 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)));

KurtosisAccumulator!(double, KurtosisAlgo.threePass, Summation.naive) w = 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;
T scaledSumOfQuarts;
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;
Summator!(T, summation) centeredSumOfQuarts;
void put(Range)(Range r) if (isIterable!Range);
void put()(T x);
void put()(KurtosisAccumulator!(T, kurtosisAlgo, summation) v);
@property size_t count();
@property F mean(F = T)();
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

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];

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);

Examples:

test basic functionality

import mir.math.common: approxEqual;
import mir.ndslice.slice: sliced;

EntropyAccumulator!(double, Summation.pairwise) x;
x.put([0.1, 0.2, 0.3].sliced);
assert(x.entropy.approxEqual(-0.913338));
x.put(0.4);
assert(x.entropy.approxEqual(-1.279854));

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)

If T is a floating point type, this is an alias to the unqualified type.

If T is not a floating point type, this will alias a double type if T is summable and implicitly convertible to a floating point type.

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:

Summation

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];

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];

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];

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 a type for which isComplex!F is true.

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 (that are either implicitly convertible to floating point or if isComplex is true)

import mir.ndslice.slice: sliced;
import std.complex: Complex;
import std.math.operations: isClose;

auto x = [Complex!double(1, 2), Complex!double(2, 3), Complex!double(3, 4), Complex!double(4, 5)].sliced;
assert(x.rawMoment!2.isClose(Complex!double(-24, 80)/ 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 a type for which isComplex!F is true.

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];

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 (that are either implicitly convertible to floating point or if isComplex is true)

import mir.ndslice.slice: sliced;
import std.complex: Complex;
import std.math.operations: isClose;

auto x = [Complex!double(1, 2), Complex!double(2, 3), Complex!double(3, 4), Complex!double(4, 5)].sliced;
assert(x.centralMoment!2.isClose(Complex!double(0, 10) / 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, then the result will have a double type if F is implicitly convertible to a floating point 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

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];

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];

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

Library Reference

mir.stat.descriptive.univariate

Discontinuous sample quantile

Continuous sample quantile