参考: 安川章, 「できる人が使っている最小二乗法の一発フィット」,『Interface(インターフェース) 2017年 08 月号』CQ出版, p.143
差の 2 乗や偏微分の計算などはせずに参考記事に従って擬似逆行列を使って計算してみる。
何かを測定して次のデータ対が得られたとする。このデータが N 次函数で近似できるものと假定する。
x = {1,2,3,4,5,6,7,8,9,10}
y = {3,5,4,2,6,7,9,7,6,4}
このデータを次のように行列で表現する (N 次 → 0 次の順に並べたのは、多項式を組み立てる polyEval() 函数が高次 → 低次の順に引数を取るからという理由に過ぎない)。
a = [ [x[0]^N, x[0]^(N-1), ......, x[0]^0], [x[1]^N, x[1]^(N-1), ......, x[1]^0], . . . [x[k]^N, x[k]^(N-1), ......, x[k]^0] ] b = [ [y[0]], [y[1]], ......, [y[k]] ] coeffs = [ [coeffs[N]], [coeffs[N-1]], ......, [coeffs[0]] ] a * coeffs = b
未知数 coeffs = b * (a^-1) が近似多項式の各次数の係数である。理窟はわからないがこれで最小二乗法と同じ結果になるとの由。
一般にデータのほうが未知数よりも圧倒的に多いため上の行列 a は正方行列にはならず逆行列 (a^-1) が求まらない。そこで擬似逆行列を利用する。擬似逆行列は次式で求まる。
a の擬似逆行列 = ((a の転置行列) * a)^-1 * (a の転置行列)
.tns (least_square.tns - Google ドライブ)
Define least_square(xlist, ylist, order)= Func :Local a, a_trans, a_pinv, coeffs, poly_func :a:= seq(seq(xlist[m]^(n), n, order, 0, −1), m, 1, dim(xlist)) © 上の行列 a を組み立てる。 :a_trans:= a@t © 行列 a の転置行列を求める。 :a_pinv:= (a_trans * a)^(−1) * a_trans © 行列 a の擬似逆行列を求める。 :b:= (list▶mat(ylist))@t © 上の行列 b を組み立てる。 :coeffs:= a_pinv * b © 擬似逆行列を使って多項式の各次数の係数を求める。 :coeffs:= mat▶list(coeffs) :poly_func:= polyEval(coeffs, x) © 近似多項式を組み立てて、 :Return approx(poly_func) © 返す。 :EndFunc
実行結果: