Skip to content

Commit 41c6f44

Browse files
insel-mazinsel-maz
committed
Add multilinear regression
1 parent c064e80 commit 41c6f44

File tree

8 files changed

+664
-19
lines changed

8 files changed

+664
-19
lines changed

‎src/LinearAlgebra/MatrixFactory.php‎

Lines changed: 24 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -5,6 +5,7 @@
55
use MathPHP\Exception;
66
use MathPHP\Number\Complex;
77
use MathPHP\Number\ObjectArithmetic;
8+
use MathPHP\Polynomials\MonomialExponentGenerator;
89

910
/**
1011
* Matrix factory to create matrices of all types.
@@ -524,7 +525,7 @@ public static function hilbert(int $n): NumericMatrix
524525
/**
525526
* Create the Vandermonde Matrix from a simple array.
526527
*
527-
* @param array<int|float> $M (α₁, α₂, α₃ ⋯ αm)
528+
* @param array<int|float>|array<array<int|float>> $M (α₁, α₂, α₃ ⋯ αm)
528529
* @param int $n
529530
*
530531
* @return NumericMatrix
@@ -537,9 +538,28 @@ public static function hilbert(int $n): NumericMatrix
537538
public static function vandermonde(array $M, int $n): NumericMatrix
538539
{
539540
$A = [];
540-
foreach ($M as $row => $α) {
541-
for ($i = 0; $i < $n; $i++) {
542-
$A[$row][$i] = $α ** $i;
541+
if (!empty($M)) {
542+
// Create at least a one-column matrix.
543+
$M = \array_map(function ($value) {
544+
return \is_array($value) ? $value : [$value];
545+
}, $M);
546+
547+
$dimension = \count(\reset($M));
548+
$degree = $n - 1;
549+
$exponentTuples = MonomialExponentGenerator::all($dimension, $degree, true);
550+
551+
foreach ($M as $row) {
552+
$values = [];
553+
foreach ($exponentTuples as $exponents) {
554+
$value = 1; // start as int
555+
\reset($row);
556+
foreach ($exponents as $exponent) {
557+
$value *= \current($row) ** $exponent;
558+
\next($row);
559+
}
560+
$values[] = $value;
561+
}
562+
$A[] = $values;
543563
}
544564
}
545565

‎src/Polynomials/MonomialExponentGenerator.php‎

Lines changed: 101 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,101 @@
1+
<?php
2+
3+
namespace MathPHP\Polynomials;
4+
5+
use Generator;
6+
use InvalidArgumentException;
7+
8+
final class MonomialExponentGenerator
9+
{
10+
/**
11+
* Returns all exponent tuples with total degree <= $degree.
12+
*
13+
* @param int $dimension d >= 1
14+
* @param int $degree p >= 0
15+
* @param bool $reverse
16+
* @return list<list<int>>
17+
*/
18+
public static function all(int $dimension, int $degree, bool $reverse): array
19+
{
20+
$gen = self::iterate($dimension, $degree, $reverse);
21+
return \iterator_to_array($gen, false);
22+
}
23+
24+
/**
25+
* Generator over all exponent tuples with total degree <= $degree.
26+
* Uses generators to keep memory usage low for large d/p.
27+
*
28+
* @param int $dimension
29+
* @param int $degree
30+
* @param bool $reverse
31+
* @return Generator<int, list<int>> yields int[] (length = $dimension)
32+
*/
33+
public static function iterate(int $dimension, int $degree, bool $reverse): Generator
34+
{
35+
if ($dimension < 1) {
36+
throw new InvalidArgumentException("dimension must be >= 1.");
37+
}
38+
if ($degree < 0) {
39+
throw new InvalidArgumentException("degree must be >= 0.");
40+
}
41+
42+
$current = \array_fill(0, $dimension, 0);
43+
44+
if ($reverse) {
45+
// Degrees 0..p; within each degree use lexicographic order
46+
for ($g = 0; $g <= $degree; $g++) {
47+
yield from self::recursiveDistributeRevLex($dimension, $g, 0, $current);
48+
}
49+
} else {
50+
// Degrees 0..p; within each degree use lexicographic order
51+
for ($g = 0; $g <= $degree; $g++) {
52+
yield from self::recursiveDistributeLex($dimension, $g, 0, $current);
53+
}
54+
}
55+
}
56+
57+
/**
58+
* Recursive helper: distributes `remaining` units across positions in lexicographic order.
59+
*
60+
* @param int $dimension
61+
* @param int $remaining
62+
* @param int $pos
63+
* @param list<int> $current Variable reference for performance.
64+
* @return Generator<int, list<int>>
65+
*/
66+
private static function recursiveDistributeLex(int $dimension, int $remaining, int $pos, array &$current): Generator
67+
{
68+
if ($pos === $dimension - 1) {
69+
$current[$pos] = $remaining;
70+
yield $current;
71+
return;
72+
}
73+
for ($e = 0; $e <= $remaining; $e++) {
74+
$current[$pos] = $e;
75+
yield from self::recursiveDistributeLex($dimension, $remaining - $e, $pos + 1, $current);
76+
}
77+
}
78+
79+
/**
80+
* Recursive helper: distributes `remaining` units across positions in reverse-lex order.
81+
*
82+
* @param int $dimension
83+
* @param int $remaining
84+
* @param int $pos
85+
* @param list<int> $current Variable reference for performance.
86+
* @return Generator<int, list<int>>
87+
*/
88+
private static function recursiveDistributeRevLex(int $dimension, int $remaining, int $pos, array &$current): Generator
89+
{
90+
if ($pos === $dimension - 1) {
91+
$current[$pos] = $remaining;
92+
yield $current;
93+
return;
94+
}
95+
// reverse-lex: prioritize larger exponents at earlier positions
96+
for ($e = $remaining; $e >= 0; $e--) {
97+
$current[$pos] = $e;
98+
yield from self::recursiveDistributeRevLex($dimension, $remaining - $e, $pos + 1, $current);
99+
}
100+
}
101+
}

‎src/Statistics/Regression/Methods/LeastSquares.php‎

Lines changed: 3 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -25,10 +25,10 @@ trait LeastSquares
2525
private $reg_ys;
2626

2727
/**
28-
* Regression xs
28+
* Regression xs or xss
2929
* Since the actual xs may be translated for regression, we need to keep these
3030
* handy for regression statistics.
31-
* @var array<float>
31+
* @var array<float>|array<array<float>>
3232
*/
3333
private $reg_xs;
3434

@@ -96,7 +96,7 @@ trait LeastSquares
9696
* (x)² - x²
9797
*
9898
* @param array<float> $ys y values
99-
* @param array<float> $xs x values
99+
* @param array<float>|array<array<float>> $xs x values
100100
* @param int $order The order of the polynomial. 1 = linear, 2 = x², etc
101101
* @param int $fit_constant '1' if we are fitting a constant to the regression.
102102
*

‎src/Statistics/Regression/Models/MultilinearModel.php‎

Lines changed: 79 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,79 @@
1+
<?php
2+
3+
namespace MathPHP\Statistics\Regression\Models;
4+
5+
trait MultilinearModel
6+
{
7+
private static $subscripts = ['₀', '₁', '₂', '₃', '₄', '₅', '₆', '₇', '₈', '₉'];
8+
9+
/**
10+
* @param int $n
11+
*
12+
* @return string
13+
*/
14+
public static function getSubscript(int $n): string
15+
{
16+
return \implode(
17+
'',
18+
\array_map(
19+
function ($c) {
20+
return self::$subscripts[$c];
21+
},
22+
\str_split((string)$n)
23+
)
24+
);
25+
}
26+
27+
/**
28+
* Evaluate the model given all the model parameters
29+
* y = x₁β₁ + x₂β₂ + … + xₖβₖ + ε
30+
*
31+
* @param list<float> $vector
32+
* @param non-empty-list<float> $params
33+
*
34+
* @return float y evaluated
35+
*/
36+
public static function evaluateModel(array $vector, array $params): float
37+
{
38+
$y = 0;
39+
for ($i = 1, $len = \count($params); $i < $len; $i++) {
40+
$y += $vector[$i - 1] * $params[$i];
41+
}
42+
$y += $params[0];
43+
return $y;
44+
}
45+
46+
/**
47+
* Get regression parameters (coefficients)
48+
*
49+
* @param non-empty-list<float> $params
50+
*
51+
* @return array<string, float>
52+
*/
53+
public function getModelParameters(array $params): array
54+
{
55+
$result = [];
56+
for ($i = 1, $len = \count($params); $i < $len; $i++) {
57+
$result['β' . self::getSubscript($i)] = $params[$i];
58+
}
59+
$result['ε'] = $params[0];
60+
return $result;
61+
}
62+
63+
/**
64+
* Get regression equation
65+
*
66+
* @param non-empty-list<float> $params
67+
*
68+
* @return string
69+
*/
70+
public function getModelEquation(array $params): string
71+
{
72+
$result = 'y = ';
73+
for ($i = 1, $len = \count($params); $i < $len; $i++) {
74+
$result .= \sprintf('%fx%s + ', $params[$i], self::getSubscript($i));
75+
}
76+
$result .= \sprintf('%f', $params[0]);
77+
return $result;
78+
}
79+
}

‎src/Statistics/Regression/Multilinear.php‎

Lines changed: 52 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,52 @@
1+
<?php
2+
3+
namespace MathPHP\Statistics\Regression;
4+
5+
use MathPHP\Exception;
6+
use MathPHP\Exception\BadDataException;
7+
8+
class Multilinear extends ParametricRegression
9+
{
10+
use Methods\LeastSquares;
11+
use Models\MultilinearModel;
12+
13+
/**
14+
* Calculates the regression parameters.
15+
*
16+
* @throws Exception\BadDataException
17+
* @throws Exception\IncorrectTypeException
18+
* @throws Exception\MatrixException
19+
* @throws Exception\MathException
20+
*/
21+
public function calculate(): void
22+
{
23+
$this->parameters = $this->leastSquares($this->ys, $this->xss)->getColumn(0);
24+
}
25+
26+
/**
27+
* Evaluate the regression equation at x.
28+
*
29+
* @param float $x
30+
*
31+
* @return float
32+
*
33+
* @throws BadDataException
34+
*/
35+
public function evaluate(float $x): float
36+
{
37+
throw new Exception\BadDataException('Multilinear regression does not support evaluate(x)');
38+
}
39+
40+
/**
41+
* Evaluate the regression equation at x vector.
42+
* Uses the instance model's evaluateModel method.
43+
*
44+
* @param non-empty-list<float> $vector
45+
*
46+
* @return float
47+
*/
48+
public function evaluateVector(array $vector): float
49+
{
50+
return $this->evaluateModel($vector, $this->parameters);
51+
}
52+
}

0 commit comments

Comments
 (0)