juliaで最小二乗法を解く

ファイルの読み込みに少し悪戦苦闘しつつも最小二乗法が完成した。

次のようなファイルデータから最小二乗法で線形回帰を行う。

x	y
5.6	30
5.8	26
6.0	33
6.2	31
6.4	33
6.4	35
6.4	37
6.6	36
6.8	33

として、一番いいとを求める。

番目のデータを, としたときの次のを最小化するとを求める。
\[
J=\frac{1}{2}\sum_{i=1}^n(y_i - (ax_i+b))^2
\]
それぞれ偏微分すると次のようになる。
\[
\begin{cases}
\frac{\partial J}{\partial a} = \sum_{i=1}^n(y_i - (ax_i+b))(-x_i) \\
\frac{\partial J}{\partial b} = \sum_{i=1}^n(y_i - (ax_i+b))(-1)
\end{cases}
\]
これらをとおいて整理すると次のようになるから、左から逆行列を掛ければ解ける。
\[
\left(
\begin{array}{cc}
\sum_{i=1}^nx_i^2 & \sum_{i=i}^nx_i \\
\sum_{i=1}^nx_i & n \\
\end{array}
\right)
\left(
\begin{array}{cc}
a \\
b \\
\end{array}
\right) =
\left(
\begin{array}{cc}
\sum_{i=1}^nx_iy_i \\
\sum_{i=1}^ny_i \\
\end{array}
\right)
\]

io = open("data.txt", "r")
n = countlines(io) #サンプル数
data = zeros(Float64, n, 2) # n*2行列

seekstart(io)
for i = 1:n
    line = split(readline(io))
    data[i, 1] = parse(Float64, line[1])
    data[i, 2] = parse(Float64, line[2])
end

x = data[:,1]
y = data[:,2]
X = sum(x) # Σxi
Y = sum(y) # Σyi
X2 = x'x # Σ(xi^2)
XY = x'y # Σxiyi

A = [X2 X; X n]
B = [XY; Y]

print(A\B) # A^-1 * B = [6.03383, -5.01128]

本当に1からjuliaを触っているのでしょうもないことで時間が取られてる印象。
次は多変量。