Adding kde() method for Kernel Density Estimation

Author: roberto-butti

CHANGELOG.md

@@ -1,5 +1,8 @@
 # Changelog
 
+## 1.2.3 - WIP
+- Adding `kde()` method for Kernel Density Estimation — returns a closure that estimates PDF or CDF from sample data, supporting 9 kernel functions with aliases
+
 ## 1.2.2 - 2026-02-21
 - Adding `method` parameter to `quantiles()` supporting `'exclusive'` (default) and `'inclusive'` interpolation methods
 - Adding `medianGrouped()` method for estimating the median of grouped/binned continuous data using interpolation

README.md

@@ -82,6 +82,7 @@ The various mathematical statistics are listed below:
 | `correlation()` | Pearson’s or Spearman’s rank correlation coefficient for two inputs |
 | `covariance()` | the sample covariance of two inputs |
 | `linearRegression()` | return the slope and intercept of simple linear regression parameters estimated using ordinary least squares (supports `proportional: true` for regression through the origin) |
+| `kde()` | kernel density estimation — returns a closure that estimates the probability density (or CDF) at any point |
 
 #### Stat::mean( array $data )
 Return the sample arithmetic mean of the array _$data_.
@@ -400,6 +401,34 @@ list($slope, $intercept) = Stat::linearRegression(
 // $intercept = 0.0
 ```
 
+#### Stat::kde ( array $data , float $h , string $kernel = 'normal' , bool $cumulative = false )
+Create a continuous probability density function (or cumulative distribution function) from discrete sample data using Kernel Density Estimation.
+Returns a `Closure` that can be called with any point to estimate the density (or CDF value).
+
+Supported kernels: `normal` (alias `gauss`), `logistic`, `sigmoid`, `rectangular` (alias `uniform`), `triangular`, `parabolic` (alias `epanechnikov`), `quartic` (alias `biweight`), `triweight`, `cosine`.
+
+```php
+$data = [-2.1, -1.3, -0.4, 1.9, 5.1, 6.2];
+$f = Stat::kde($data, h: 1.5);
+$f(2.5);
+// estimated density at x = 2.5
+```
+
+Using a different kernel:
+
+```php
+$f = Stat::kde($data, h: 1.5, kernel: 'triangular');
+$f(2.5);
+```
+
+Cumulative distribution function:
+
+```php
+$F = Stat::kde($data, h: 1.5, cumulative: true);
+$F(2.5);
+// estimated CDF at x = 2.5 (probability that a value is <= 2.5)
+```
+
 ### Freq class
 With *Statistics* package you can calculate frequency table.
 A frequency table lists the frequency of various outcomes in a sample.

src/Stat.php

@@ -715,6 +715,171 @@ private static function ranks(array $data): array
         return $ranks;
     }
 
+    /**
+     * Create a continuous probability density function or cumulative distribution
+     * function from discrete sample data using Kernel Density Estimation.
+     *
+     * Returns a Closure that estimates the density (or CDF) at any given point.
+     *
+     * @param  array<int|float>  $data  sample data
+     * @param  float  $h  bandwidth (smoothing parameter), must be > 0
+     * @param  string  $kernel  kernel name (normal, logistic, sigmoid, rectangular, triangular, parabolic, quartic, triweight, cosine) or alias
+     * @param  bool  $cumulative  if true, return CDF estimator; otherwise PDF estimator
+     * @return \Closure  a callable that takes a float and returns the estimated density or CDF value
+     *
+     * @throws InvalidDataInputException if data is empty, bandwidth <= 0, or kernel is invalid
+     */
+    public static function kde(
+        array $data,
+        float $h,
+        string $kernel = 'normal',
+        bool $cumulative = false,
+    ): \Closure {
+        if ($data === []) {
+            throw new InvalidDataInputException("The data must not be empty.");
+        }
+        if ($h <= 0) {
+            throw new InvalidDataInputException("Bandwidth h must be positive.");
+        }
+
+        $aliases = [
+            'gauss' => 'normal',
+            'uniform' => 'rectangular',
+            'epanechnikov' => 'parabolic',
+            'biweight' => 'quartic',
+        ];
+        $kernel = strtolower($kernel);
+        if (isset($aliases[$kernel])) {
+            $kernel = $aliases[$kernel];
+        }
+
+        $sqrt2pi = sqrt(2.0 * M_PI);
+
+        // Standard normal CDF using Abramowitz & Stegun approximation (7.1.26)
+        $normalCdf = static function (float $t) use ($sqrt2pi): float {
+            $negative = $t < 0;
+            $t = abs($t);
+            $b1 = 0.319381530;
+            $b2 = -0.356563782;
+            $b3 = 1.781477937;
+            $b4 = -1.821255978;
+            $b5 = 1.330274429;
+            $p = 0.2316419;
+            $k = 1.0 / (1.0 + $p * $t);
+            $pdf = exp(-$t * $t / 2.0) / $sqrt2pi;
+            $cdf = 1.0 - $pdf * $k * ($b1 + $k * ($b2 + $k * ($b3 + $k * ($b4 + $k * $b5))));
+
+            return $negative ? 1.0 - $cdf : $cdf;
+        };
+
+        $kernels = [
+            'normal' => [
+                'pdf' => static fn(float $t): float => exp(-$t * $t / 2.0) / $sqrt2pi,
+                'cdf' => $normalCdf,
+                'support' => null,
+            ],
+            'logistic' => [
+                'pdf' => static fn(float $t): float => 0.5 / (1.0 + cosh($t)),
+                'cdf' => static fn(float $t): float => 1.0 / (1.0 + exp(-$t)),
+                'support' => null,
+            ],
+            'sigmoid' => [
+                'pdf' => static fn(float $t): float => (1.0 / M_PI) / cosh($t),
+                'cdf' => static fn(float $t): float => (2.0 / M_PI) * atan(exp($t)),
+                'support' => null,
+            ],
+            'rectangular' => [
+                'pdf' => static fn(float $t): float => 0.5,
+                'cdf' => static fn(float $t): float => 0.5 * $t + 0.5,
+                'support' => 1.0,
+            ],
+            'triangular' => [
+                'pdf' => static fn(float $t): float => 1.0 - abs($t),
+                'cdf' => static fn(float $t): float => $t >= 0
+                    ? 1.0 - (1.0 - $t) * (1.0 - $t) / 2.0
+                    : (1.0 + $t) * (1.0 + $t) / 2.0,
+                'support' => 1.0,
+            ],
+            'parabolic' => [
+                'pdf' => static fn(float $t): float => 0.75 * (1.0 - $t * $t),
+                'cdf' => static fn(float $t): float => -0.25 * $t * $t * $t + 0.75 * $t + 0.5,
+                'support' => 1.0,
+            ],
+            'quartic' => [
+                'pdf' => static fn(float $t): float => (15.0 / 16.0) * (1.0 - $t * $t) ** 2,
+                'cdf' => static fn(float $t): float => (15.0 * $t - 10.0 * $t ** 3 + 3.0 * $t ** 5) / 16.0 + 0.5,
+                'support' => 1.0,
+            ],
+            'triweight' => [
+                'pdf' => static fn(float $t): float => (35.0 / 32.0) * (1.0 - $t * $t) ** 3,
+                'cdf' => static fn(float $t): float => (35.0 * $t - 35.0 * $t ** 3 + 21.0 * $t ** 5 - 5.0 * $t ** 7) / 32.0 + 0.5,
+                'support' => 1.0,
+            ],
+            'cosine' => [
+                'pdf' => static fn(float $t): float => (M_PI / 4.0) * cos(M_PI * $t / 2.0),
+                'cdf' => static fn(float $t): float => 0.5 * sin(M_PI * $t / 2.0) + 0.5,
+                'support' => 1.0,
+            ],
+        ];
+
+        if (! isset($kernels[$kernel])) {
+            $valid = implode(', ', array_merge(array_keys($kernels), array_keys($aliases)));
+            throw new InvalidDataInputException(
+                "Unknown kernel '{$kernel}'. Valid kernels: {$valid}.",
+            );
+        }
+
+        $kernelDef = $kernels[$kernel];
+        $support = $kernelDef['support'];
+        $fn = $cumulative ? $kernelDef['cdf'] : $kernelDef['pdf'];
+
+        $sorted = $data;
+        sort($sorted);
+        $n = count($sorted);
+
+        if ($cumulative) {
+            return static function (float $x) use ($sorted, $n, $h, $fn, $support): float {
+                $sum = 0.0;
+                if ($support !== null) {
+                    $lo = self::bisectLeft($sorted, $x - $h * $support);
+                    $hi = self::bisectRight($sorted, $x + $h * $support);
+                    for ($i = $lo; $i < $hi; $i++) {
+                        $t = ($x - $sorted[$i]) / $h;
+                        $sum += $fn($t);
+                    }
+                    // Points entirely to the left contribute 1.0 each
+                    $sum += $lo;
+                } else {
+                    for ($i = 0; $i < $n; $i++) {
+                        $t = ($x - $sorted[$i]) / $h;
+                        $sum += $fn($t);
+                    }
+                }
+
+                return $sum / $n;
+            };
+        }
+
+        return static function (float $x) use ($sorted, $n, $h, $fn, $support): float {
+            $sum = 0.0;
+            if ($support !== null) {
+                $lo = self::bisectLeft($sorted, $x - $h * $support);
+                $hi = self::bisectRight($sorted, $x + $h * $support);
+                for ($i = $lo; $i < $hi; $i++) {
+                    $t = ($x - $sorted[$i]) / $h;
+                    $sum += $fn($t);
+                }
+            } else {
+                for ($i = 0; $i < $n; $i++) {
+                    $t = ($x - $sorted[$i]) / $h;
+                    $sum += $fn($t);
+                }
+            }
+
+            return $sum / ($n * $h);
+        };
+    }
+
     /**
      * @param  array<int|float>  $x
      * @param  array<int|float>  $y

tests/StatTest.php

@@ -575,4 +575,114 @@ public function test_proportional_linear_regression_with_all_zeros_x(): void
         $this->expectException(InvalidDataInputException::class);
         Stat::linearRegression([0, 0, 0, 0, 0], [1, 2, 3, 4, 5], proportional: true);
     }
+
+    public function test_kde_normal(): void
+    {
+        $data = [-2.1, -1.3, -0.4, 1.9, 5.1, 6.2];
+        $f = Stat::kde($data, 1.5);
+        $this->assertIsCallable($f);
+
+        $density = $f(2.5);
+        $this->assertIsFloat($density);
+        $this->assertGreaterThan(0, $density);
+
+        // Verify against manually computed value
+        // f(2.5) = (1/(6*1.5)) * sum of K((2.5 - xi)/1.5)
+        $n = count($data);
+        $h = 1.5;
+        $sum = 0.0;
+        $sqrt2pi = sqrt(2.0 * M_PI);
+        foreach ($data as $xi) {
+            $t = (2.5 - $xi) / $h;
+            $sum += exp(-$t * $t / 2.0) / $sqrt2pi;
+        }
+        $expected = $sum / ($n * $h);
+        $this->assertEqualsWithDelta($expected, $density, 1e-10);
+    }
+
+    public function test_kde_all_kernels(): void
+    {
+        $data = [1.0, 2.0, 3.0, 4.0, 5.0];
+        $kernels = [
+            'normal', 'logistic', 'sigmoid',
+            'rectangular', 'triangular', 'parabolic',
+            'quartic', 'triweight', 'cosine',
+        ];
+
+        foreach ($kernels as $kernel) {
+            $f = Stat::kde($data, 1.0, $kernel);
+            $this->assertIsCallable($f, "Kernel '$kernel' should return a callable");
+            $density = $f(3.0);
+            $this->assertGreaterThanOrEqual(0, $density, "Kernel '$kernel' density should be >= 0");
+        }
+    }
+
+    public function test_kde_cumulative(): void
+    {
+        $data = [1.0, 2.0, 3.0, 4.0, 5.0];
+        $F = Stat::kde($data, 1.0, 'normal', cumulative: true);
+        $this->assertIsCallable($F);
+
+        // CDF should be monotonically non-decreasing
+        $prev = $F(-100.0);
+        foreach ([-10.0, 0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 10.0, 100.0] as $x) {
+            $current = $F($x);
+            $this->assertGreaterThanOrEqual($prev, $current, "CDF should be non-decreasing at x=$x");
+            $prev = $current;
+        }
+
+        // Approaches 0 far left and 1 far right
+        $this->assertEqualsWithDelta(0.0, $F(-100.0), 0.01);
+        $this->assertEqualsWithDelta(1.0, $F(100.0), 0.01);
+    }
+
+    public function test_kde_aliases(): void
+    {
+        $data = [1.0, 2.0, 3.0, 4.0, 5.0];
+        $x = 3.0;
+
+        // gauss == normal
+        $f1 = Stat::kde($data, 1.0, 'gauss');
+        $f2 = Stat::kde($data, 1.0, 'normal');
+        $this->assertEqualsWithDelta($f1($x), $f2($x), 1e-15);
+
+        // uniform == rectangular
+        $f1 = Stat::kde($data, 1.0, 'uniform');
+        $f2 = Stat::kde($data, 1.0, 'rectangular');
+        $this->assertEqualsWithDelta($f1($x), $f2($x), 1e-15);
+
+        // epanechnikov == parabolic
+        $f1 = Stat::kde($data, 1.0, 'epanechnikov');
+        $f2 = Stat::kde($data, 1.0, 'parabolic');
+        $this->assertEqualsWithDelta($f1($x), $f2($x), 1e-15);
+
+        // biweight == quartic
+        $f1 = Stat::kde($data, 1.0, 'biweight');
+        $f2 = Stat::kde($data, 1.0, 'quartic');
+        $this->assertEqualsWithDelta($f1($x), $f2($x), 1e-15);
+    }
+
+    public function test_kde_empty_data(): void
+    {
+        $this->expectException(InvalidDataInputException::class);
+        Stat::kde([], 1.0);
+    }
+
+    public function test_kde_invalid_bandwidth(): void
+    {
+        $this->expectException(InvalidDataInputException::class);
+        Stat::kde([1.0, 2.0], 0.0);
+    }
+
+    public function test_kde_invalid_bandwidth_negative(): void
+    {
+        $this->expectException(InvalidDataInputException::class);
+        Stat::kde([1.0, 2.0], -1.0);
+    }
+
+    public function test_kde_invalid_kernel(): void
+    {
+        $this->expectException(InvalidDataInputException::class);
+        Stat::kde([1.0, 2.0], 1.0, 'invalid_kernel');
+    }
 }