Introduce Welford's online algorithm for variance (#116)

Author: sanmai

.github/workflows/main.yml

@@ -113,6 +113,10 @@ jobs:
         run: |
           make yamllint
 
+      - name: Non-code tests
+        run: |
+          php vendor/bin/phpunit --group=documentation
+
       - name: Static Analysis
         run: |
           make ci-analyze --keep-going

README.md

@@ -125,6 +125,8 @@ All entry points always return an instance of the pipeline.
 | `max()`     | Finds the highest value. | `max` |
 | `min()`     | Finds the lowest value. | `min` |
 | `toArray()` | Returns an array with all values. Eagerly executed. | `dict`, `ToDictionary` |
+| `runningVariance()` | Computes online statistics: sample mean, sample variance, standard deviation. | [Welford's method](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm) |
+| `finalVariance()` | Computes final statistics for the sequence. |   |
 | `__construct()` | Can be provided with an optional initial iterator. Used in the `take()` function from above.  |     |
 
 Pipeline is an iterator and can be used as any other iterable. 
@@ -394,6 +396,66 @@ foreach ($pipeline as $value) {
 
 This allows to skip type checks for return values if one has no results to return: instead of `false` or `null` it is safe to return an unprimed pipeline.
 
+## `$pipeline->runningVariance()`
+
+Computes online statistics for the sequence: counts, sample  mean, sample variance, standard deviation. You can access these numbers on the fly with methods such as `getCount()`, `getMean()`, `getVariance()`, `getStandardDeviation()`.
+
+This method also accepts an optional cast callback that should return `float|null`: `null` values are discarded. Therefore you can have several running variances computing numbers for different parts of the data.
+
+```php
+$pipeline->runningVariance($varianceForShippedOrders, static function (order $order): ?float {
+    if (!$order->isShipped()) {
+        // This order will be excluded from the computation.
+        return null;
+    }
+
+    return $order->getTotal();
+});
+
+$pipeline->runningVariance($varianceForPaidOrders, static function (order $order): ?float {
+    if ($order->isUnpaid()) {
+        // This order will be excluded from the computation.
+        return null;
+    }
+
+    return $order->getProjectedTotal();
+});
+```
+
+As you process orders, you will be able to access `$varianceForShippedOrders->getMean()` and `$varianceForPaidOrders->getMean()`.
+
+This computation uses [Welford's online algorithm](https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm), therefore it can handle very large numbers of data points.
+
+## `$pipeline->finalVariance()`
+
+A convenience method to computes the final statistics for the sequence. Accepts an optional cast method, else assumes the sequence countains valid numbers.
+
+```php
+// Fibonacci numbers generator
+$fibonacci = map(function () {
+    yield 0;
+
+    $prev = 0;
+    $current = 1;
+
+    while (true) {
+        yield $current;
+        $next = $prev + $current;
+        $prev = $current;
+        $current = $next;
+    }
+});
+
+// Statistics for the second hundred Fibonacci numbers.
+$variance = $fibonacci->slice(101, 100)->finalVariance();
+
+$variance->getStandardDeviation();
+// float(3.5101061922557E+40)
+
+$variance->getCount();
+// int(100)
+```
+
 # Contributions
 
 Contributions to documentation and test cases are welcome. Bug reports are welcome too. 

example.php

@@ -66,6 +66,9 @@
     yield $i * $j;
 });
 
+// Collect online statistics
+$pipeline->runningVariance($variance);
+
 // one way to filter
 $pipeline->map(function ($i) {
     if ($i > 50) {
@@ -87,6 +90,13 @@
 var_dump($value);
 // int(104)
 
+var_dump($variance->getCount());
+// int(12)
+var_dump($variance->getMean());
+// float(22)
+var_dump($variance->getStandardDeviation());
+// float(30.3794188704967)
+
 $sum = take([1, 2, 3])->reduce();
 var_dump($sum);
 // int(6)
@@ -137,3 +147,27 @@
 //     [1] =>
 //     int(3)
 // }
+
+// Fibonacci numbers generator
+$fibonacci = map(function () {
+    yield 0;
+
+    $prev = 0;
+    $current = 1;
+
+    while (true) {
+        yield $current;
+        $next = $prev + $current;
+        $prev = $current;
+        $current = $next;
+    }
+});
+
+// Statistics for the second hundred Fibonacci numbers
+$variance = $fibonacci->slice(101, 100)->finalVariance();
+
+var_dump($variance->getStandardDeviation());
+// float(3.5101061922557E+40)
+
+var_dump($variance->getCount());
+// int(100)

infection.json.dist

@@ -11,6 +11,7 @@
         ],
         "@default": true,
         "IdenticalEqual": false,
-        "NotIdenticalNotEqual": false
+        "NotIdenticalNotEqual": false,
+        "DecrementInteger": false
     }
 }

phpunit.xml.dist

@@ -18,10 +18,9 @@
         </include>
     </source>
 
-    <coverage>
-        <report>
-          <clover outputFile="build/logs/clover.xml"/>
-          <text outputFile="php://stdout"/>
-        </report>
-    </coverage>
+    <groups>
+        <exclude>
+            <group>documentation</group>
+        </exclude>
+    </groups>
 </phpunit>

psalm.xml.dist

@@ -20,6 +20,7 @@
                 <file name="example.php" />
             </errorLevel>
         </ForbiddenCode>
+
         <MissingClosureParamType errorLevel="suppress" />
         <MissingClosureReturnType errorLevel="suppress" />
 

src/Helper/RunningVariance.php

@@ -0,0 +1,141 @@
+<?php
+/**
+ * Copyright 2017, 2018 Alexey Kopytko <[email protected]>
+ *
+ * Licensed under the Apache License, Version 2.0 (the "License");
+ * you may not use this file except in compliance with the License.
+ * You may obtain a copy of the License at
+ *
+ * http://www.apache.org/licenses/LICENSE-2.0
+ *
+ * Unless required by applicable law or agreed to in writing, software
+ * distributed under the License is distributed on an "AS IS" BASIS,
+ * WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+ * See the License for the specific language governing permissions and
+ * limitations under the License.
+ */
+
+declare(strict_types=1);
+
+namespace Pipeline\Helper;
+
+use function sqrt;
+
+/**
+ * Computes statistics (such as standard deviation) in real time.
+ *
+ * @see https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Welford's_online_algorithm
+ *
+ * @final
+ */
+class RunningVariance
+{
+    /**
+     * The number of observed values.
+     *
+     * @var int<0, max>
+     */
+    private int $count = 0;
+
+    /** First moment: the mean value. */
+    private float $mean = 0.0;
+
+    /** Second moment: the aggregated squared distance from the mean. */
+    private float $m2 = 0.0;
+
+    public function __construct(self ...$spiesToMerge)
+    {
+        foreach ($spiesToMerge as $spy) {
+            $this->merge($spy);
+        }
+    }
+
+    public function observe(float $value): float
+    {
+        ++$this->count;
+
+        $delta = $value - $this->mean;
+
+        $this->mean += $delta / $this->count;
+        $this->m2 += $delta * ($value - $this->mean);
+
+        return $value;
+    }
+
+    /**
+     * The number of observed values.
+     */
+    public function getCount(): int
+    {
+        return $this->count;
+    }
+
+    /**
+     * Get the mean value.
+     */
+    public function getMean(): float
+    {
+        if (0 === $this->count) {
+            // For no values the variance is undefined.
+            return NAN;
+        }
+
+        return $this->mean;
+    }
+
+    /**
+     * Get the variance.
+     */
+    public function getVariance(): float
+    {
+        if (0 === $this->count) {
+            // For no values the variance is undefined.
+            return NAN;
+        }
+
+        if (1 === $this->count) {
+            // Avoiding division by zero: variance for one value is zero.
+            return 0.0;
+        }
+
+        // https://en.wikipedia.org/wiki/Bessel%27s_correction
+        return $this->m2 / ($this->count - 1);
+    }
+
+    /**
+     * Compute the standard deviation.
+     */
+    public function getStandardDeviation(): float
+    {
+        return sqrt($this->getVariance());
+    }
+
+    /**
+     * Merge another instance into this instance.
+     *
+     * @see https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
+     */
+    private function merge(self $other): void
+    {
+        // Shortcut a no-op
+        if (0 === $other->count) {
+            return;
+        }
+
+        // Avoid division by zero by copying values
+        if (0 === $this->count) {
+            $this->count = $other->count;
+            $this->mean = $other->mean;
+            $this->m2 = $other->m2;
+
+            return;
+        }
+
+        $count = $this->count + $other->count;
+        $delta = $other->mean - $this->mean;
+
+        $this->mean = ($this->count * $this->mean) / $count + ($other->count * $other->mean) / $count;
+        $this->m2 = $this->m2 + $other->m2 + ($delta ** 2 * $this->count * $other->count / $count);
+        $this->count = $count;
+    }
+}

src/Standard.php

@@ -987,4 +987,73 @@ private static function flipKeysAndValues(iterable $previous): iterable
             yield $value => $key;
         }
     }
+
+    private function feedRunningVariance(Helper\RunningVariance $variance, ?callable $castFunc): self
+    {
+        if (null === $castFunc) {
+            $castFunc = 'floatval';
+        }
+
+        return $this->cast(static function ($value) use ($variance, $castFunc) {
+            $float = $castFunc($value);
+
+            if (null !== $float) {
+                $variance->observe($float);
+            }
+
+            // Returning the original value here
+            return $value;
+        });
+    }
+
+    /**
+     * Feeds in an instance of RunningVariance.
+     *
+     * @param ?Helper\RunningVariance &$variance the instance of RunningVariance; initialized unless provided
+     * @param ?callable               $castFunc  the cast callback, returning ?float; null values are not counted
+     *
+     * @param-out Helper\RunningVariance $variance
+     *
+     * @return $this
+     */
+    public function runningVariance(
+        ?Helper\RunningVariance &$variance,
+        ?callable $castFunc = null
+    ): self {
+        $variance ??= new Helper\RunningVariance();
+
+        $this->feedRunningVariance($variance, $castFunc);
+
+        return $this;
+    }
+
+    /**
+     * Computes final statistics for the sequence.
+     *
+     * @param ?callable               $castFunc the cast callback, returning ?float; null values are not counted
+     * @param ?Helper\RunningVariance $variance the optional instance of RunningVariance
+     */
+    public function finalVariance(
+        ?callable $castFunc = null,
+        ?Helper\RunningVariance $variance = null
+    ): Helper\RunningVariance {
+        $variance ??= new Helper\RunningVariance();
+
+        if ($this->empty()) {
+            // No-op: an empty array.
+            return $variance;
+        }
+
+        $this->feedRunningVariance($variance, $castFunc);
+
+        if (is_array($this->pipeline)) {
+            // We are done!
+            return $variance;
+        }
+
+        // Consume every available item
+        $_ = iterator_count($this->pipeline);
+
+        return $variance;
+    }
 }

tests/DocumentationTest.php

@@ -35,6 +35,8 @@
 /**
  * Test documentation has sections for all public methods.
  *
+ * @group documentation
+ *
  * @coversNothing
  *
  * @internal