Add multilinear regression

src/LinearAlgebra/MatrixFactory.php

@@ -5,6 +5,7 @@
 use MathPHP\Exception;
 use MathPHP\Number\Complex;
 use MathPHP\Number\ObjectArithmetic;
+use MathPHP\Polynomials\MonomialExponentGenerator;
 
 /**
  * Matrix factory to create matrices of all types.
@@ -524,7 +525,7 @@ public static function hilbert(int $n): NumericMatrix
     /**
      * Create the Vandermonde Matrix from a simple array.
      *
-     * @param array<int|float> $M (α₁, α₂, α₃ ⋯ αm)
+     * @param array<int|float>|array<array<int|float>> $M (α₁, α₂, α₃ ⋯ αm)
      * @param int   $n
      *
      * @return NumericMatrix
@@ -537,9 +538,28 @@ public static function hilbert(int $n): NumericMatrix
     public static function vandermonde(array $M, int $n): NumericMatrix
     {
         $A = [];
-        foreach ($M as $row => $α) {
-            for ($i = 0; $i < $n; $i++) {
-                $A[$row][$i] = $α ** $i;
+        if (!empty($M)) {
+            // Create at least a one-column matrix.
+            $M = \array_map(function ($value) {
+                return \is_array($value) ? $value : [$value];
+            }, $M);
+
+            $dimension = \count(\reset($M));
+            $degree = $n - 1;
+            $exponentTuples = MonomialExponentGenerator::all($dimension, $degree, true);
+
+            foreach ($M as $row) {
+                $values = [];
+                foreach ($exponentTuples as $exponents) {
+                    $value = 1; // start as int
+                    \reset($row);
+                    foreach ($exponents as $exponent) {
+                        $value *= \current($row) ** $exponent;
+                        \next($row);
+                    }
+                    $values[] = $value;
+                }
+                $A[] = $values;
             }
         }
 

src/Polynomials/MonomialExponentGenerator.php

@@ -0,0 +1,115 @@
+<?php
+
+namespace MathPHP\Polynomials;
+
+use Generator;
+use InvalidArgumentException;
+use MathPHP\Exception;
+use MathPHP\Probability\Combinatorics;
+
+final class MonomialExponentGenerator
+{
+    /**
+     * @param int $dimension
+     * @param int $degree
+     * @return int
+     * @throws Exception\OutOfBoundsException
+     */
+    public static function getNumberOfTerms(int $dimension, int $degree): int
+    {
+        return (int)(Combinatorics::factorial($dimension + $degree) /
+            (Combinatorics::factorial($degree) * Combinatorics::factorial($dimension)));
+    }
+
+    /**
+     * Returns all exponent tuples with total degree <= $degree.
+     *
+     * @param int $dimension d >= 1
+     * @param int $degree p >= 0
+     * @param bool $reverse
+     * @return list<list<int>>
+     */
+    public static function all(int $dimension, int $degree, bool $reverse): array
+    {
+        $gen = self::iterate($dimension, $degree, $reverse);
+        return \iterator_to_array($gen, false);
+    }
+
+    /**
+     * Generator over all exponent tuples with total degree <= $degree.
+     * Uses generators to keep memory usage low for large d/p.
+     *
+     * @param int $dimension
+     * @param int $degree
+     * @param bool $reverse
+     * @return Generator<int, list<int>> yields int[] (length = $dimension)
+     */
+    public static function iterate(int $dimension, int $degree, bool $reverse): Generator
+    {
+        if ($dimension < 1) {
+            throw new InvalidArgumentException("dimension must be >= 1.");
+        }
+        if ($degree < 0) {
+            throw new InvalidArgumentException("degree must be >= 0.");
+        }
+
+        $current = \array_fill(0, $dimension, 0);
+
+        if ($reverse) {
+            // Degrees 0..p; within each degree use lexicographic order
+            for ($g = 0; $g <= $degree; $g++) {
+                yield from self::recursiveDistributeRevLex($dimension, $g, 0, $current);
+            }
+        } else {
+            // Degrees 0..p; within each degree use lexicographic order
+            for ($g = 0; $g <= $degree; $g++) {
+                yield from self::recursiveDistributeLex($dimension, $g, 0, $current);
+            }
+        }
+    }
+
+    /**
+     * Recursive helper: distributes `remaining` units across positions in lexicographic order.
+     *
+     * @param int $dimension
+     * @param int $remaining
+     * @param int $pos
+     * @param list<int> $current Variable reference for performance.
+     * @return Generator<int, list<int>>
+     */
+    private static function recursiveDistributeLex(int $dimension, int $remaining, int $pos, array &$current): Generator
+    {
+        if ($pos === $dimension - 1) {
+            $current[$pos] = $remaining;
+            yield $current;
+            return;
+        }
+        for ($e = 0; $e <= $remaining; $e++) {
+            $current[$pos] = $e;
+            yield from self::recursiveDistributeLex($dimension, $remaining - $e, $pos + 1, $current);
+        }
+    }
+
+    /**
+     * Recursive helper: distributes `remaining` units across positions in reverse-lex order.
+     *
+     * @param int $dimension
+     * @param int $remaining
+     * @param int $pos
+     * @param list<int> $current Variable reference for performance.
+     * @return Generator<int, list<int>>
+     */
+    private static function recursiveDistributeRevLex(int $dimension, int $remaining, int $pos, array &$current): Generator
+    {
+        if ($pos === $dimension - 1) {
+            $current[$pos] = $remaining;
+            yield $current;
+            return;
+        }
+        // reverse-lex: prioritize larger exponents at earlier positions
+        for ($e = $remaining; $e >= 0; $e--) {
+            $current[$pos] = $e;
+            yield from self::recursiveDistributeRevLex($dimension, $remaining - $e, $pos + 1, $current);
+        }
+    }
+}

src/Statistics/Regression/HanesWoolf.php

@@ -12,12 +12,25 @@
  *       K + x
  *
  * The equation is linearized and fit using Least Squares
+ *
+ * @phpstan-import-type SimpleLinearResultModel from Methods\LeastSquares
+ * @phpstan-import-type PolynomialResultModel from Methods\LeastSquares
  */
 class HanesWoolf extends ParametricRegression
 {
+    /** @use Methods\LeastSquares<SimpleLinearResultModel> */
     use Methods\LeastSquares;
     use Models\MichaelisMenten;
 
+    /**
+     * @param list<float> $array
+     * @return SimpleLinearResultModel
+     */
+    protected function createResultModel(array $array): array
+    {
+        return $this->createSimpleLinearResultModel($array);
+    }
+
     /**
      * Calculate the regression parameters by least squares on linearized data
      * x / y = x / V + K / V

src/Statistics/Regression/Linear.php

@@ -24,12 +24,25 @@
  *
  *     _    _
  * b = y - mx
+ *
+ * @phpstan-import-type SimpleLinearResultModel from Methods\LeastSquares
+ * @phpstan-import-type PolynomialResultModel from Methods\LeastSquares
  */
 class Linear extends ParametricRegression
 {
+    /** @use Methods\LeastSquares<SimpleLinearResultModel> */
     use Methods\LeastSquares;
     use Models\LinearModel;
 
+    /**
+     * @param list<float> $array
+     * @return SimpleLinearResultModel
+     */
+    protected function createResultModel(array $array): array
+    {
+        return $this->createSimpleLinearResultModel($array);
+    }
+
     /**
      * Calculates the regression parameters.
      *

src/Statistics/Regression/LinearThroughPoint.php

@@ -23,9 +23,13 @@
  *         ∑(x-v)²
  *
  * b = w - m * v
+ *
+ * @phpstan-import-type SimpleLinearResultModel from Methods\LeastSquares
+ * @phpstan-import-type PolynomialResultModel from Methods\LeastSquares
  */
 class LinearThroughPoint extends ParametricRegression
 {
+    /** @use Methods\LeastSquares<SimpleLinearResultModel> */
     use Methods\LeastSquares;
     use Models\LinearModel;
 
@@ -53,6 +57,15 @@ public function __construct(array $points, array $force = [0,0])
         parent::__construct($points);
     }
 
+    /**
+     * @param list<float> $array
+     * @return SimpleLinearResultModel
+     */
+    protected function createResultModel(array $array): array
+    {
+        return $this->createSimpleLinearResultModel($array);
+    }
+
     /**
      * Calculates the regression parameters.
      *

src/Statistics/Regression/LineweaverBurk.php

@@ -12,12 +12,25 @@
  *       K + x
  *
  * The equation is linearized and fit using Least Squares
+ *
+ * @phpstan-import-type SimpleLinearResultModel from Methods\LeastSquares
+ * @phpstan-import-type PolynomialResultModel from Methods\LeastSquares
  */
 class LineweaverBurk extends ParametricRegression
 {
     use Models\MichaelisMenten;
+    /** @use Methods\LeastSquares<SimpleLinearResultModel> */
     use Methods\LeastSquares;
 
+    /**
+     * @param list<float> $array
+     * @return SimpleLinearResultModel
+     */
+    protected function createResultModel(array $array): array
+    {
+        return $this->createSimpleLinearResultModel($array);
+    }
+
     /**
      * Calculate the regression parameters by least squares on linearized data
      * y⁻¹ = K * V⁻¹ * x⁻¹ + V⁻¹