【入門】多項式回帰分析(Julia)【数値計算】

【入門】多項式回帰分析(Julia)【数値計算】 数値計算
【入門】多項式回帰分析(Julia)【数値計算】

MATLAB、Python、Scilab、Julia比較ページはこちら
https://www.simulationroom999.com/blog/comparison-of-matlab-python-scilab/

はじめに

の、
MATLAB,Python,Scilab,Julia比較 第2章 その66【多項式回帰分析⑤】

正規方程式を用いた、多項式回帰分析について。
今回は、Juliaで演算してみる。

正規方程式、各パラメータ、推定対象の多項式再掲

まずは正規方程式、多項式回帰分析に於ける各パラメータ、推定対象の多項式の再掲。

正規方程式

\(
x=(A^TA)^{-1}A^Tb
\)

多項式回帰分析に於ける各パラメータ

\(
A=
\begin{bmatrix}
x_1^2 & x_1 & 1\\
x_2^2 & x_2 & 1\\
\vdots & \vdots & \vdots\\
x_n^2 & x_n & 1\\
\end{bmatrix},
\vec{x}=
\begin{bmatrix}
\alpha\\
\beta\\
\gamma
\end{bmatrix},
\vec{b}=
\begin{bmatrix}
y_1\\
y_2\\
\vdots\\
y_n
\end{bmatrix}
\)

推定対象の多項式

\(
z=4x^2-5y+2
\)

これをJuliaで解いてみる。

Juliaコード

Juliaコードは以下になる。

using PyPlot

n = 100

x = rand(1, n)
y = 4*x.^2 -5*x.+2+rand(1, n).-0.5
A=[x'.^2   x'   ones(length(x),1)]
b=y'
X=(A'*A)^-1 *A'*b
print(X)

plot(x, y, "+b")
xp=range(0, 1, length=100);
yp=range(0, 1, length=100);

plot( xp, X[1]*xp.^2+X[2]*xp.+X[3], "r",linewidth=3)

処理結果

処理結果は以下。

正規方程式で多項式回帰分析(Julia)、Figure 1
正規方程式で多項式回帰分析(Julia)、Figure 1
[3.9143559582970897; -5.107059369648499; 2.0972036620536327;;]

考察

Juliaの挙動としてはOK。
これもコード自体はMATLABからのほぼコピペ。
ただし、linspaceをrangeにしたり、plotのオプションがちょっと違ったりで微調整はした。
(なまじ似てるが故にケアレスミスもしやすい・・・。)

まとめ

  • 正規方程式による多項式回帰分析をJuliaで実施。
  • 誤差はあるものの目的の係数の算出はできている。
  • コード自体はMATLABコードのほぼコピペ。
    • 等差数列、plotのオプション周りの合わせこみはした。

バックナンバーはこちら。

コメント

タイトルとURLをコピーしました