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,101 @@
+<?php
+
+namespace MathPHP\Polynomials;
+
+use Generator;
+use InvalidArgumentException;
+
+final class MonomialExponentGenerator
+{
+    /**
+     * 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/Methods/LeastSquares.php

@@ -25,10 +25,10 @@ trait LeastSquares
     private $reg_ys;
 
     /**
-     * Regression xs
+     * Regression xs or xss
      * Since the actual xs may be translated for regression, we need to keep these
      * handy for regression statistics.
-     * @var array<float>
+     * @var array<float>|array<array<float>>
      */
     private $reg_xs;
 
@@ -96,7 +96,7 @@ trait LeastSquares
      *       (x)² - x²
      *
      * @param  array<float> $ys y values
-     * @param  array<float> $xs x values
+     * @param  array<float>|array<array<float>> $xs x values
      * @param  int   $order The order of the polynomial. 1 = linear, 2 = x², etc
      * @param  int   $fit_constant '1' if we are fitting a constant to the regression.
      *

src/Statistics/Regression/Models/MultilinearModel.php

@@ -0,0 +1,81 @@
+<?php
+
+namespace MathPHP\Statistics\Regression\Models;
+
+trait MultilinearModel
+{
+    /**
+     * @var list<string>
+     */
+    private static $subscripts = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'];
+
+    /**
+     * @param int $n
+     *
+     * @return string
+     */
+    public static function getSubscript(int $n): string
+    {
+        return \implode(
+            '',
+            \array_map(
+                function ($c) {
+                    return self::$subscripts[$c];
+                },
+                \str_split((string)$n)
+            )
+        );
+    }
+
+    /**
+     * Evaluate the model given all the model parameters
+     * y = β₀ + x₁β₁ + x₂β₂ + … + xₖβₖ
+     *
+     * @param list<float> $vector
+     * @param non-empty-list<float> $params
+     *
+     * @return float y evaluated
+     */
+    public static function evaluateModel(array $vector, array $params): float
+    {
+        $y = $params[0];
+        for ($i = 1, $len = \count($params); $i < $len; $i++) {
+            $y += $vector[$i - 1] * $params[$i];
+        }
+        return $y;
+    }
+
+    /**
+     * Get regression parameters (coefficients).
+     *
+     * Use array_values() on the result to get a list<float> of coefficients.
+     *
+     * @param array<int, float> $params
+     *
+     * @return array<string, float> [β₀ + β₁ + β₂ + … + βₖ]
+     */
+    public function getModelParameters(array $params): array
+    {
+        $result = [];
+        for ($i = 0, $len = \count($params); $i < $len; $i++) {
+            $result['β' . self::getSubscript($i)] = $params[$i];
+        }
+        return $result;
+    }
+
+    /**
+     * Get regression equation
+     *
+     * @param non-empty-list<float> $params
+     *
+     * @return string
+     */
+    public function getModelEquation(array $params): string
+    {
+        $result = \sprintf('y = %f', $params[0]);
+        for ($i = 1, $len = \count($params); $i < $len; $i++) {
+            $result .= \sprintf(' + %fx%s', $params[$i], self::getSubscript($i));
+        }
+        return $result;
+    }
+}

src/Statistics/Regression/Multilinear.php

@@ -0,0 +1,57 @@
+<?php
+
+namespace MathPHP\Statistics\Regression;
+
+use MathPHP\Exception;
+use MathPHP\Exception\BadDataException;
+
+class Multilinear extends ParametricRegression
+{
+    use Methods\LeastSquares;
+    use Models\MultilinearModel;
+
+    /**
+     * Calculates the regression parameters.
+     *
+     * @throws Exception\BadDataException
+     * @throws Exception\IncorrectTypeException
+     * @throws Exception\MatrixException
+     * @throws Exception\MathException
+     */
+    public function calculate(): void
+    {
+        $this->parameters = $this->leastSquares($this->ys, $this->xss)->getColumn(0);
+    }
+
+    /**
+     * Evaluate the regression equation at x.
+     *
+     * @param float $x
+     *
+     * @return float
+     *
+     * @throws BadDataException
+     */
+    public function evaluate(float $x): float
+    {
+        throw new Exception\BadDataException('Multilinear regression does not support evaluate(x)');
+    }
+
+    /**
+     * Evaluate the regression equation at x vector.
+     * Uses the instance model's evaluateModel method.
+     *
+     * @param non-empty-list<float> $vector
+     *
+     * @return float
+     *
+     * @throws BadDataException
+     */
+    public function evaluateVector(array $vector): float
+    {
+        if (empty($this->parameters)) {
+            throw new Exception\BadDataException('Regression parameters are not calculated');
+        }
+        return $this->evaluateModel($vector, $this->parameters);
+    }
+}