diff --git a/lib/pure/stats.nim b/lib/pure/stats.nim new file mode 100644 index 0000000000..a78a2b5bf5 --- /dev/null +++ b/lib/pure/stats.nim @@ -0,0 +1,269 @@ +# +# +# Nim's Runtime Library +# (c) Copyright 2015 Andreas Rumpf +# +# See the file "copying.txt", included in this +# distribution, for details about the copyright. +# +## Statistical analysis framework for performing +## basic statistical analysis of data. +## The data is analysed in a single pass, when a data value +## is pushed to the ``RunningStat`` or ``RunningRegress`` objects +## +## ``RunningStat`` calculates for a single data set +## - n (data count) +## - min (smallest value) +## - max (largest value) +## - mean +## - variance +## - varianceS (sample var) +## - standardDeviation +## - standardDeviationS (sample stddev) +## - skewness (the third statistical moment) +## - kurtosis (the fourth statistical moment) +## +## ``RunningRegress`` calculates for two sets of data +## - n +## - slope +## - intercept +## - correlation +## +from math import FloatClass, sqrt, pow, round + +{.push debugger:off .} # the user does not want to trace a part + # of the standard library! +{.push checks:off, line_dir:off, stack_trace:off.} + +type + RunningStat* = object ## an accumulator for statistical data + n*: int ## number of pushed data + min*, max*: float ## self-explaining + mom1, mom2, mom3, mom4: float ## statistical moments, mom1 is mean + + + RunningRegress* = object ## an accumulator for regression calculations + n*: int ## number of pushed data + x_stats*: RunningStat ## stats for first set of data + y_stats*: RunningStat ## stats for second set of data + s_xy: float ## acculated data for combined xy + +{.deprecated: [TFloatClass: FloatClass, TRunningStat: RunningStat].} + +proc clear(s: var RunningStat) = + ## reset `s` + s.n = 0 + s.min = 0.0 + s.max = 0.0 + s.mom1 = 0.0 + s.mom2 = 0.0 + s.mom3 = 0.0 + s.mom4 = 0.0 + +proc push*(s: var RunningStat, x: float) = + ## pushes a value `x` for processing + inc(s.n) + # See Knuth TAOCP vol 2, 3rd edition, page 232 + if s.min > x: s.min = x + if s.max < x: s.max = x + let n = toFloat(s.n) + let delta = x - s.mom1 + let delta_n = delta / toFloat(s.n) + let delta_n2 = delta_n * delta_n + let term1 = delta * delta_n * toFloat(s.n - 1) + s.mom4 += term1 * delta_n2 * (n*n - 3*n + 3) + + 6*delta_n2*s.mom2 - 4*delta_n*s.mom3 + s.mom3 += term1 * delta_n * (n - 2) - 3*delta_n*s.mom2 + s.mom2 += term1 + s.mom1 += delta_n + +proc push*(s: var RunningStat, x: int) = + ## pushes a value `x` for processing. + ## + ## `x` is simply converted to ``float`` + ## and the other push operation is called. + s.push(toFloat(x)) + +proc push*(s: var RunningStat, x: openarray[float|int]) = + ## pushes all values of `x` for processing. + ## + ## Int values of `x` are simply converted to ``float`` and + ## the other push operation is called. + for val in x: + s.push(val) + +proc mean*(s: RunningStat): float = result = s.mom1 + ## computes the current mean of `s` + +proc variance*(s: RunningStat): float = + ## computes the current population variance of `s` + result = s.mom2 / toFloat(s.n) + +proc varianceS*(s: RunningStat): float = + ## computes the current sample variance of `s` + if s.n > 1: result = s.mom2 / toFloat(s.n - 1) + +proc standardDeviation*(s: RunningStat): float = + ## computes the current population standard deviation of `s` + result = sqrt(variance(s)) + +proc standardDeviationS*(s: RunningStat): float = + ## computes the current sample standard deviation of `s` + result = sqrt(varianceS(s)) + +proc skewness*(s: RunningStat): float = + ## computes the current population skewness of `s` + result = sqrt(toFloat(s.n)) * s.mom3 / pow(s.mom2, 1.5) + +proc skewnessS*(s: RunningStat): float = + ## computes the current sample skewness of `s` + let s2 = skewness(s) + result = sqrt(toFloat(s.n*(s.n-1)))*s2 / toFloat(s.n-2) + +proc kurtosis*(s: RunningStat): float = + ## computes the current population kurtosis of `s` + result = toFloat(s.n) * s.mom4 / (s.mom2 * s.mom2) - 3.0 + +proc kurtosisS*(s: RunningStat): float = + ## computes the current sample kurtosis of `s` + result = toFloat(s.n-1) / toFloat((s.n-2)*(s.n-3)) * + (toFloat(s.n+1)*kurtosis(s) + 6) + +proc `+`*(a, b: RunningStat): RunningStat = + ## combine two RunningStats. + ## + ## Useful if performing parallel analysis of data series + ## and need to re-combine parallel result sets + result.clear() + result.n = a.n + b.n + + let delta = b.mom1 - a.mom1 + let delta2 = delta*delta + let delta3 = delta*delta2 + let delta4 = delta2*delta2 + let n = toFloat(result.n) + + result.mom1 = (a.n.float*a.mom1 + b.n.float*b.mom1) / n + result.mom2 = a.mom2 + b.mom2 + delta2 * a.n.float * b.n.float / n + result.mom3 = a.mom3 + b.mom3 + + delta3 * a.n.float * b.n.float * (a.n.float - b.n.float)/(n*n); + result.mom3 += 3.0*delta * (a.n.float*b.mom2 - b.n.float*a.mom2) / n + result.mom4 = a.mom4 + b.mom4 + + delta4*a.n.float*b.n.float * toFloat(a.n*a.n - a.n*b.n + b.n*b.n) / + (n*n*n) + result.mom4 += 6.0*delta2 * (a.n.float*a.n.float*b.mom2 + b.n.float*b.n.float*a.mom2) / + (n*n) + + 4.0*delta*(a.n.float*b.mom3 - b.n.float*a.mom3) / n + result.max = max(a.max, b.max) + result.min = max(a.min, b.min) + +proc `+=`*(a: var RunningStat, b: RunningStat) {.inline.} = + ## add a second RunningStats `b` to `a` + a = a + b + +# ---------------------- Running Regression ----------------------------- + +proc clear*(r: var RunningRegress) = + ## reset `r` + r.x_stats.clear() + r.y_stats.clear() + r.s_xy = 0.0 + r.n = 0 + +proc push*(r: var RunningRegress, x, y: float) = + ## pushes two values `x` and `y` for processing + r.s_xy += (r.x_stats.mean() -x)*(r.y_stats.mean() - y)*toFloat(r.n)/toFloat(r.n+1) + + r.x_stats.push(x) + r.y_stats.push(y) + inc(r.n) + +proc push*(r: var RunningRegress, x, y: int) {.inline.} = + ## pushes two values `x` and `y` for processing. + ## + ## `x` and `y` are converted to ``float`` + ## and the other push operation is called. + r.push(toFloat(x), toFloat(y)) + +proc push*(r: var RunningRegress, x, y: openarray[float|int]) = + ## pushes two sets of values `x` and `y` for processing. + assert(x.len == y.len) + for i in 0..