|
| 1 | +/++ |
| 2 | +This module contains base statistical algorithms. |
| 3 | +
|
| 4 | +Note that used specialized summing algorithms execute more primitive operations |
| 5 | +than vanilla summation. Therefore, if in certain cases maximum speed is required |
| 6 | +at expense of precision, one can use $(SUBREF, sum, Summation.fast). |
| 7 | +
|
| 8 | +License: $(LINK2 http://boost.org/LICENSE_1_0.txt, Boost License 1.0). |
| 9 | +
|
| 10 | +Authors: Ilya Yaroshenko |
| 11 | +
|
| 12 | +Copyright: 2019 Symmetry Investments Group and Kaleidic Associates Advisory Limited. |
| 13 | +
|
| 14 | +Macros: |
| 15 | +SUBREF = $(REF_ALTTEXT $(TT $2), $2, mir, ndslice, $1)$(NBSP) |
| 16 | +T2=$(TR $(TDNW $(LREF $1)) $(TD $+)) |
| 17 | +T4=$(TR $(TDNW $(LREF $1)) $(TD $2) $(TD $3) $(TD $4)) |
| 18 | ++/ |
| 19 | +module mir.math.stat; |
| 20 | + |
| 21 | +import core.lifetime: move; |
| 22 | +import mir.math.common: optmath; |
| 23 | +import mir.math.sum; |
| 24 | +import mir.primitives; |
| 25 | +import std.range.primitives: isInputRange; |
| 26 | +import std.traits: isArray, isFloatingPoint; |
| 27 | + |
| 28 | +/++ |
| 29 | +Computes the average of `r`, which must be a finite iterable. |
| 30 | +
|
| 31 | +Returns: |
| 32 | + The average of all the elements in the range r. |
| 33 | ++/ |
| 34 | +template mean(Summation summation = Summation.appropriate) |
| 35 | +{ |
| 36 | + /// |
| 37 | + @safe @optmath sumType!Range |
| 38 | + mean(Range)(Range r) |
| 39 | + if (hasLength!Range |
| 40 | + || summation == Summation.appropriate |
| 41 | + || summation == Summation.fast |
| 42 | + || summation == Summation.naive) |
| 43 | + { |
| 44 | + static if (hasLength!Range) |
| 45 | + { |
| 46 | + auto n = r.length; |
| 47 | + return sum!summation(r.move) / cast(sumType!Range) n; |
| 48 | + } |
| 49 | + else |
| 50 | + { |
| 51 | + auto s = cast(typeof(return)) 0; |
| 52 | + size_t length; |
| 53 | + foreach (e; r) |
| 54 | + { |
| 55 | + length++; |
| 56 | + s += e; |
| 57 | + } |
| 58 | + return s / cast(sumType!Range) length; |
| 59 | + } |
| 60 | + } |
| 61 | +} |
| 62 | + |
| 63 | +///ditto |
| 64 | +template mean(string summation) |
| 65 | +{ |
| 66 | + mixin("alias mean = .mean!(Summation." ~ summation ~ ");"); |
| 67 | +} |
| 68 | + |
| 69 | +/// |
| 70 | +version(mir_test) @safe pure nothrow unittest |
| 71 | +{ |
| 72 | + assert(mean([1.0, 2, 3]) == 2); |
| 73 | +} |
| 74 | + |
| 75 | +/++ |
| 76 | +A linear regression model with a single explanatory variable. |
| 77 | ++/ |
| 78 | +template simpleLinearRegression(Summation summation = Summation.kbn) |
| 79 | +{ |
| 80 | + import mir.ndslice.slice; |
| 81 | + |
| 82 | + /++ |
| 83 | + Params: |
| 84 | + x = `x[i]` points |
| 85 | + y = `f(x[i])` values |
| 86 | + Returns: |
| 87 | + The pair of shift and slope of the linear curve. |
| 88 | + +/ |
| 89 | + @optmath |
| 90 | + sumType!YRange[2] |
| 91 | + simpleLinearRegression(XRange, YRange)(XRange x, YRange y) @safe |
| 92 | + if (isInputRange!XRange && isInputRange!YRange && !(isArray!XRange && isArray!YRange) && isFloatingPoint!(sumType!YRange)) |
| 93 | + in { |
| 94 | + static if (hasLength!XRange && hasLength!YRange) |
| 95 | + assert(x.length == y.length); |
| 96 | + } |
| 97 | + body { |
| 98 | + alias X = typeof(sumType!XRange.init * sumType!XRange.init); |
| 99 | + alias Y = sumType!YRange; |
| 100 | + enum summationX = !__traits(isIntegral, X) ? summation : Summation.naive; |
| 101 | + Summator!(X, summationX) xms = 0; |
| 102 | + Summator!(Y, summation) yms = 0; |
| 103 | + Summator!(X, summationX) xxms = 0; |
| 104 | + Summator!(Y, summation) xyms = 0; |
| 105 | + |
| 106 | + static if (hasLength!XRange) |
| 107 | + sizediff_t n = x.length; |
| 108 | + else |
| 109 | + sizediff_t n = 0; |
| 110 | + |
| 111 | + while (!x.empty) |
| 112 | + { |
| 113 | + static if (!(hasLength!XRange && hasLength!YRange)) |
| 114 | + assert(!y.empty); |
| 115 | + |
| 116 | + static if (!hasLength!XRange) |
| 117 | + n++; |
| 118 | + |
| 119 | + auto xi = x.front; |
| 120 | + auto yi = y.front; |
| 121 | + xms.put(xi); |
| 122 | + yms.put(yi); |
| 123 | + xxms.put(xi * xi); |
| 124 | + xyms.put(xi * yi); |
| 125 | + |
| 126 | + y.popFront; |
| 127 | + x.popFront; |
| 128 | + } |
| 129 | + |
| 130 | + static if (!(hasLength!XRange && hasLength!YRange)) |
| 131 | + assert(y.empty); |
| 132 | + |
| 133 | + auto xm = xms.sum; |
| 134 | + auto ym = yms.sum; |
| 135 | + auto xxm = xxms.sum; |
| 136 | + auto xym = xyms.sum; |
| 137 | + |
| 138 | + auto slope = (xym * n - xm * ym) / (xxm * n - xm * xm); |
| 139 | + |
| 140 | + return [(ym - slope * xm) / n, slope]; |
| 141 | + } |
| 142 | + |
| 143 | + /// ditto |
| 144 | + @optmath |
| 145 | + typeof(X.init * Y.init)[2] |
| 146 | + simpleLinearRegression(X, Y)(scope const X[] x, scope const Y[] y) @safe |
| 147 | + { |
| 148 | + return .simpleLinearRegression!summation(x.sliced, y.sliced); |
| 149 | + } |
| 150 | +} |
| 151 | + |
| 152 | +/// ditto |
| 153 | +template simpleLinearRegression(string summation) |
| 154 | +{ |
| 155 | + mixin("alias simpleLinearRegression = .simpleLinearRegression!(Summation." ~ summation ~ ");"); |
| 156 | +} |
| 157 | + |
| 158 | +/// |
| 159 | +version(mir_test) @safe pure nothrow @nogc unittest |
| 160 | +{ |
| 161 | + import mir.math.common: approxEqual; |
| 162 | + static immutable x = [0, 1, 2, 3]; |
| 163 | + static immutable y = [-1, 0.2, 0.9, 2.1]; |
| 164 | + auto params = x.simpleLinearRegression(y); |
| 165 | + assert(params[0].approxEqual(-0.95)); // shift |
| 166 | + assert(params[1].approxEqual(1)); // slope |
| 167 | +} |
0 commit comments