Add mir.math.stat

doc/Makefile

@@ -37,7 +37,7 @@ PACKAGE_mir_combinatorics = package
 PACKAGE_mir_container = binaryheap
 PACKAGE_mir_graph = tarjan package
 PACKAGE_mir_interpolate = package constant linear spline utility polynomial
-PACKAGE_mir_math = sum numeric
+PACKAGE_mir_math = sum numeric stat
 PACKAGE_mir_math_func = expdigamma
 PACKAGE_mir_ndslice_connect = cpython
 PACKAGE_mir_rc = package ptr array context

index.d

@@ -18,7 +18,8 @@ $(BOOKTABLE ,
  $(LEADINGROW Math)
  $(TR $(TDNW $(MREF mir,numeric)) $(TD Basic numeric optimisations))
  $(TR $(TDNW $(MREF mir,math,numeric)) $(TD Simple numeric algorithms))
- $(TR $(TDNW $(MREF mir,math,sum)★) $(TD Various precise summation algorithms))
+ $(TR $(TDNW $(MREF mir,math,sum)) $(TD Various precise summation algorithms))
+ $(TR $(TDNW $(MREF mir,math,stat)) $(TD Basic API for statistics))
  $(TR $(TDNW $(MREF mir,math,constant)) $(TD Math constants))
  $(TR $(TDNW $(MREF mir,polynomial)) $(TD Polynomial ref-counted structure))
  $(LEADINGROW Reference counting)

meson.build

@@ -32,6 +32,7 @@ mir_algorithm_src = [
  'source/mir/interpolate/utility.d',
  'source/mir/math/func/expdigamma.d',
  'source/mir/math/numeric.d',
+ 'source/mir/math/stat.d',
  'source/mir/math/sum.d',
  'source/mir/ndslice/allocation.d',
  'source/mir/ndslice/chunks.d',
@@ -50,6 +51,7 @@ mir_algorithm_src = [
  'source/mir/ndslice/topology.d',
  'source/mir/ndslice/traits.d',
  'source/mir/numeric.d',
+ 'source/mir/polinomial.d',
  'source/mir/range.d',
  'source/mir/rc/array.d',
  'source/mir/rc/context.d',

source/mir/math/stat.d

@@ -0,0 +1,167 @@
+/++
+This module contains base statistical algorithms.
+
+Note that used specialized summing algorithms execute more primitive operations
+than vanilla summation. Therefore, if in certain cases maximum speed is required
+at expense of precision, one can use $(SUBREF, sum, Summation.fast).
+
+License: $(LINK2 http://boost.org/LICENSE_1_0.txt, Boost License 1.0).
+
+Authors: Ilya Yaroshenko
+
+Copyright: 2019 Symmetry Investments Group and Kaleidic Associates Advisory Limited.
+
+Macros:
+SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, ndslice, $1)$(NBSP)
+T2=$(TR $(TDNW $(LREF $1)) $(TD $+))
+T4=$(TR $(TDNW $(LREF $1)) $(TD $2) $(TD $3) $(TD $4))
++/
+module mir.math.stat;
+
+import core.lifetime: move;
+import mir.math.common: optmath;
+import mir.math.sum;
+import mir.primitives;
+import std.range.primitives: isInputRange;
+import std.traits: isArray, isFloatingPoint;
+
+/++
+Computes the average of `r`, which must be a finite iterable.
+
+Returns:
+ The average of all the elements in the range r.
++/
+template mean(Summation summation = Summation.appropriate)
+{
+ ///
+ @safe @optmath sumType!Range
+ mean(Range)(Range r)
+ if (hasLength!Range
+ || summation == Summation.appropriate
+ || summation == Summation.fast
+ || summation == Summation.naive)
+ {
+ static if (hasLength!Range)
+ {
+ auto n = r.length;
+ return sum!summation(r.move) / cast(sumType!Range) n;
+ }
+ else
+ {
+ auto s = cast(typeof(return)) 0;
+ size_t length;
+ foreach (e; r)
+ {
+ length++;
+ s += e;
+ }
+ return s / cast(sumType!Range) length;
+ }
+ }
+}
+
+///ditto
+template mean(string summation)
+{
+ mixin("alias mean = .mean!(Summation." ~ summation ~ ");");
+}
+
+///
+version(mir_test) @safe pure nothrow unittest
+{
+ assert(mean([1.0, 2, 3]) == 2);
+}
+
+/++
+A linear regression model with a single explanatory variable.
++/
+template simpleLinearRegression(Summation summation = Summation.kbn)
+{
+ import mir.ndslice.slice;
+
+ /++
+ Params:
+ x = `x[i]` points
+ y = `f(x[i])` values
+ Returns:
+ The pair of shift and slope of the linear curve.
+ +/
+ @optmath
+ sumType!YRange[2]
+ simpleLinearRegression(XRange, YRange)(XRange x, YRange y) @safe
+ if (isInputRange!XRange && isInputRange!YRange && !(isArray!XRange && isArray!YRange) && isFloatingPoint!(sumType!YRange))
+ in {
+ static if (hasLength!XRange && hasLength!YRange)
+ assert(x.length == y.length);
+ }
+ body {
+ alias X = typeof(sumType!XRange.init * sumType!XRange.init);
+ alias Y = sumType!YRange;
+ enum summationX = !__traits(isIntegral, X) ? summation : Summation.naive;
+ Summator!(X, summationX) xms = 0;
+ Summator!(Y, summation) yms = 0;
+ Summator!(X, summationX) xxms = 0;
+ Summator!(Y, summation) xyms = 0;
+
+ static if (hasLength!XRange)
+ sizediff_t n = x.length;
+ else
+ sizediff_t n = 0;
+
+ while (!x.empty)
+ {
+ static if (!(hasLength!XRange && hasLength!YRange))
+ assert(!y.empty);
+
+ static if (!hasLength!XRange)
+ n++;
+
+ auto xi = x.front;
+ auto yi = y.front;
+ xms.put(xi);
+ yms.put(yi);
+ xxms.put(xi * xi);
+ xyms.put(xi * yi);
+
+ y.popFront;
+ x.popFront;
+ }
+
+ static if (!(hasLength!XRange && hasLength!YRange))
+ assert(y.empty);
+
+ auto xm = xms.sum;
+ auto ym = yms.sum;
+ auto xxm = xxms.sum;
+ auto xym = xyms.sum;
+
+ auto slope = (xym * n - xm * ym) / (xxm * n - xm * xm);
+
+ return [(ym - slope * xm) / n, slope];
+ }
+
+ /// ditto
+ @optmath
+ typeof(X.init * Y.init)[2]
+ simpleLinearRegression(X, Y)(scope const X[] x, scope const Y[] y) @safe
+ {
+ return .simpleLinearRegression!summation(x.sliced, y.sliced);
+ }
+}
+
+/// ditto
+template simpleLinearRegression(string summation)
+{
+ mixin("alias simpleLinearRegression = .simpleLinearRegression!(Summation." ~ summation ~ ");");
+}
+
+///
+version(mir_test) @safe pure nothrow @nogc unittest
+{
+ import mir.math.common: approxEqual;
+ static immutable x = [0, 1, 2, 3];
+ static immutable y = [-1, 0.2, 0.9, 2.1];
+ auto params = x.simpleLinearRegression(y);
+ assert(params[0].approxEqual(-0.95)); // shift
+ assert(params[1].approxEqual(1)); // slope
+}

source/mir/math/sum.d

@@ -964,7 +964,7 @@ public:
  else
  static if (isPointer!Iterator && kind == Contiguous)
  {
- this.put(r.iterator[0 .. r.length]);
+ this.put(r.field);
  }
  else
  static if (summation == Summation.fast && N == 1)
@@ -1652,7 +1652,7 @@ unittest
  }
 }
 
-/**
+/++
 Sums elements of `r`, which must be a finite
 iterable.
 
@@ -1661,13 +1661,11 @@ value, but its type will be used if it is not specified.
 
 Note that these specialized summing algorithms execute more primitive operations
 than vanilla summation. Therefore, if in certain cases maximum speed is required
-at expense of precision, one can use $(XREF, numeric_summation, Summation.fast).
+at expense of precision, one can use $(LREF, Summation.fast).
 
 Returns:
  The sum of all the elements in the range r.
-
-See_Also: $(XREFMODULE, numeric_summation) contains detailed documentation and examples about available summation algorithms.
- */
++/
 template sum(F, Summation summation = Summation.appropriate)
  if (isFloatingPoint!F && isMutable!F)
 {
@@ -1964,7 +1962,7 @@ private T summationInitValue(T)()
  }
 }
 
-private template sumType(Range)
+package template sumType(Range)
 {
  import mir.ndslice.slice: isSlice, DeepElementType;
  static if (isSlice!Range)