最小二乗法（多項式近似）

<?php

/****************************/
/* 最小二乗法（多項式近似） */
/*      coded by Y.Suganuma */
/****************************/

	fscanf(STDIN, "%d %d", $m, $n);   // 多項式の次数とデータの数

	$x = array($n);
	$y = array($n);

	for ($i1 = 0; $i1 < $n; $i1++)   // データ
		fscanf(STDIN, "%lf %lf", $x[$i1], $y[$i1]);

	$z = least($m, $n, $x, $y);

	if ($z != NULL) {
		printf("結果\n");
		for ($i1 = 0; $i1 < $m+1; $i1++)
			printf("   %d 次の係数 %f\n", $m-$i1, $z[$i1]);
	}
	else
		printf("***error  逆行列を求めることができませんでした\n");

/******************************************/
/* 最初２乗法                             */
/*      m : 多項式の次数                  */
/*      n : データの数                    */
/*      x,y : データ                      */
/*      return : 多項式の係数（高次から） */
/*               エラーの場合はNULLを返す */
/******************************************/
function least($m, $n, $x, $y)
{
	$m++;
	$z = array($m);
	$w = array($m);
	for ($i1 = 0; $i1 < $m; $i1++)
		$w[$i1] = array($m+1);
	$A = array($n);

	for ($i1 = 0; $i1 < $n; $i1++) {
		$A[$i1]       = array($m);
		$A[$i1][$m-2] = $x[$i1];
		$A[$i1][$m-1] = 1.0;
		$x1           = $A[$i1][$m-2];
		$x2           = $x1;
		for ($i2 = $m-3; $i2 >= 0; $i2--) {
			$x2        *= $x1;
			$A[$i1][$i2]  = $x2;
		}
	}

	for ($i1 = 0; $i1 < $m; $i1++) {
		for ($i2 = 0; $i2 < $m; $i2++) {
			$w[$i1][$i2] = 0.0;
			for ($i3 = 0; $i3 < $n; $i3++)
				$w[$i1][$i2] += $A[$i3][$i1] * $A[$i3][$i2];
		}
	}

	for ($i1 = 0; $i1 < $m; $i1++) {
		$w[$i1][$m] = 0.0;
		for ($i2 = 0; $i2 < $n; $i2++)
			$w[$i1][$m] += $A[$i2][$i1] * $y[$i2];
	}

	$sw = gauss($w, $m, 1, 1.0e-10);

	if ($sw == 0) {
		for ($i1 = 0; $i1 < $m; $i1++)
			$z[$i1] = $w[$i1][$m];
	}
	else
		$z = NULL;

	return $z;
}

/*******************************************************/
/* 線形連立方程式を解く（逆行列を求める）              */
/*      w : 方程式の左辺及び右辺                       */
/*      n : 方程式の数                                 */
/*      m : 方程式の右辺の列の数                       */
/*      eps : 正則性を判定する規準                     */
/*      return : =0 : 正常                             */
/*               =1 : 逆行列が存在しない               */
/*******************************************************/
function gauss(&$w, $n, $m, $eps)
{
	$ind = 0;
	$nm  = $n + $m;

	for ($i1 = 0; $i1 < $n && $ind == 0; $i1++) {

		$y1 = .0;
		$m1 = $i1 + 1;
		$m2 = 0;

		for ($i2 = $i1; $i2 < $n; $i2++) {
			$y2 = abs($w[$i2][$i1]);
			if ($y1 < $y2) {
				$y1 = $y2;
				$m2 = $i2;
			}
		}

		if ($y1 < $eps)
			$ind = 1;

		else {

			for ($i2 = $i1; $i2 < $nm; $i2++) {
				$y1          = $w[$i1][$i2];
				$w[$i1][$i2] = $w[$m2][$i2];
				$w[$m2][$i2] = $y1;
			}

			$y1 = 1.0 / $w[$i1][$i1];

			for ($i2 = $m1; $i2 < $nm; $i2++)
				$w[$i1][$i2] *= $y1;

			for ($i2 = 0; $i2 < $n; $i2++) {
				if ($i2 != $i1) {
					for ($i3 = $m1; $i3 < $nm; $i3++)
						$w[$i2][$i3] -= $w[$i2][$i1] * $w[$i1][$i3];
				}
			}
		}
	}

	return($ind);
}

/*
-----------データ例１（コメント部分を除いて下さい）---------
1 10   // 多項式の次数とデータの数
0  2.450   // 以下，(x, y)
1  2.615
2  3.276
3  3.294
4  3.778
5  4.009
6  3.920
7  4.267
8  4.805
9  5.656

-------------------------データ例２-------------------------
2 6
10  28.2
15  47.0
20  44.4
26  32.8
32  20.8
40   0.8
*/

?>