バックナンバーはこちら。
https://www.simulationroom999.com/blog/compare-matlabpythonscilabjulia2-backnumber/
はじめに
1次関数最小二乗法の係数算出の式をJuliaを使用して実現。
登場人物
博識フクロウのフクさん
イラストACにて公開の「kino_k」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=iKciwKA9&area=1
エンジニア歴8年の太郎くん
イラストACにて公開の「しのみ」さんのイラストを使用しています。
https://www.ac-illust.com/main/profile.php?id=uCKphAW2&area=1
1次関数最小二乗法 算出式【再掲】
まずは恒例の算出式を再掲だね。
\(a,b\)を逆行列で算出
\(
\begin{bmatrix}
a \\
b
\end{bmatrix}=
\begin{bmatrix}
\sum x_i^2 && \sum x_i \\
\sum x_i && \sum 1
\end{bmatrix}^{-1}
\begin{bmatrix}
\sum x_i y_i \\
\sum y_i
\end{bmatrix}
\)
\(a,b\)を\(\sum\)で算出
\(
\begin{eqnarray}
\displaystyle a&=&\frac{n\sum x_i y_i – \sum x_i \sum y_i}{n\sum x_i^2 – (\sum x_i)^2} \\
\displaystyle b&=&\frac{-\sum x_i \sum x_i y_i + \sum x_i^2 \sum y_i}{n\sum x_i^2 – (\sum x_i)^2}
\end{eqnarray}
\)
まずは恒例の算出式を再掲だね。
まぁ今回に限ってはMATLABとほぼ一緒になると思うよ。
まぁベクトル、行列の演算自体は似てるもんねー。
Juliaコード
Juliaコードは以下になる。
using PyPlot
function LeastSquares_test()
x=[0.51, 0.76, 1.06, 1.41, 1.75, 1.9, 2.01, 2.15, 2.27, 2.4, 2.49, 2.59, 2.67, 2.76, 2.83, 2.89, 2.95, 3.01, 3.05, 3.11, 3.15, 3.19, 3.23, 3.28, 3.31, 3.34, 3.38, 3.4, 3.43, 3.46, 3.49, 3.51];
y=[10, 11, 12, 13, 14, 14.5, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, 28, 29, 30, 31, 32, 33, 34, 35, 36, 37, 38, 39, 40];
# Σで計算
print("Σで計算\n");
n = length(x);
denominator = n*sum(x.^2)-sum(x)^2;
a=(n*sum(x.*y)-sum(x)*sum(y))/denominator;
b=(-sum(x)*sum(x.*y)+sum(x.^2)*sum(y))/denominator;
print("a=",a,",b=",b,"\n");
# 行列計算
print("行列計算\n");
V_ab = inv([sum(x.^2) sum(x) ; sum(x) n])*[sum(x.*y) ; sum(y)];
print("a=",V_ab[1],",b=",V_ab[2],"\n");
xp = range(0, 4-0.01, step=0.01);
plot(x, y, marker="+", linestyle="None" );
plot(xp, a.*xp.+b );
ylim([10,41]);
xlim([0,4]);
return;
end
LeastSquares_test();結果
そして結果は以下になる。
Σで計算
a=10.133033511230964,b=-2.1616643669284774
行列計算
a=10.13303351123096,b=-2.1616643669284876これもfit関数と同じと解釈できる結果が得られたかな。
そうだね。
MATLABと比べた際のコードの差も、
linespaceがrange、plot時のオプション指定方法が若干変わるくらいか。
そうそう。
よって、今回はさほど悩まされずにすんだね。
まとめ
まとめだよ。
- 1次関数最小二乗法の係数算出の式を元にJuliaで実装。
- fit関数と同じと解釈できる結果が得られた。
- 純粋なベクトル、行列の演算に関してはMATLABとほぼ同じ書き方になる。
バックナンバーはこちら。
コメント