一元多次线性回归模型公式： y = b0 + b1* x + b2 * x^2 + b3 * x^3 + ..... + bn * x^n;

这个matlab中已经给出了方法polyfit();

下面是用php代码实现的：

public function polyfit($n, $x, $y, $poly_n, &$a)
{
for ($i = 0; $i < $n; $i ++) {
$tempx[$i] = 1;
$tempy[$i] = $y[$i];
}
for ($i = 0; $i < 2 * $poly_n + 1; $i++) {
for ($sumxx[$i] = 0, $j = 0; $j < $n; $j++)
{
$sumxx[$i] += $tempx[$j];
$tempx[$j] *= $x[$j];
}
}
for ($i = 0; $i < $poly_n + 1; $i ++) {
for ($sumxy[$i] = 0, $j = 0; $j < $n; $j ++)
{
$sumxy[$i] += $tempy[$j];
$tempy[$j] *= $x[$j];
}
}
for ($i = 0; $i < $poly_n + 1; $i ++) {
for ($j = 0; $j < $poly_n + 1; $j ++) {
$ata[$i*($poly_n + 1) + $j] = $sumxx[$i + $j];
}
}
$this->gauss_solve($poly_n + 1, $ata, $a, $sumxy);
}

public function gauss_solve($n, $A, &$x, $b)
{
for ($k = 0; $k < $n - 1; $k ++)
{
$max = abs($A[$k * $n + $k]); /*find maxmum*/
$r = $k;
for ($i = $k + 1; $i < $n - 1; $i ++) {
if ($max < abs($A[$i * $n + $i]))
{
$max = abs($A[$i * $n + $i]);
$r = $i;
}
}
if ($r != $k) {
for ($i = 0; $i < $n; $i ++) /*change array:A[k]&A[r]*/
{
$max = $A[$k * $n + $i];
$A[$k * $n + $i] = $A[$r * $n + $i];
$A[$r * $n + $i] = $max;
}
}
$max = $b[$k]; /*change array:b[k]&b[r]*/
$b[$k] = $b[$r];
$b[$r] = $max;
for ($i = $k + 1; $i < $n; $i ++)
{
for ($j = $k + 1; $j < $n; $j ++) {
$A[$i * $n + $j] -= $A[$i * $n + $k] * $A[$k * $n + $j] / $A[$k * $n + $k];
}
$b[$i] -= $A[$i * $n + $k] * $b[$k] / $A[$k * $n + $k];
}
}
for ($i = $n - 1; $i >= 0; $x[$i] /= $A[$i * $n + $i], $i --) {
for ($j = $i + 1, $x[$i] = $b[$i]; $j < $n; $j ++) {
$x[$i] -= $A[$i * $n + $j] * $x[$j];
}
}
}

 第一个函数的参数$n, $x, $y, $poly_n, &$a对应的分别是： $n: 也就是第二个参数$x的length； $x 模拟的x的值组成的数组 类似[1, 2, 3, 4, 5]; $y 模拟的y的值组成的数组； $poly_n 模拟的公式的最高次，也就是公式中的n； $a  = [] 空数组；也就是最后要求的b组成的数组；

我参考的是 http://zhidao.baidu.com/link?url=kJ5_vUgIbUhOEYafBzyqE6BDgHJIutumflEPat-_LOnwBPUcgqooMcabGVjt2WlyvPCYGkNlr49Azd77DUcUiS8c552E-55eqCsWoU3SY6e