🗃

Math＆Julia #03｜回帰モデル » 線形回帰 » 閉形式解

に公開
Julia
math
tech

はじめに

散布図にプロットされたデータを見ると、全体を貫く 1つの傾向が見えることがあります。この傾向を 1本の直線で捉えて、未知のデータの予測に活用したいとします。それでは「この傾向に最も当てはまりの良い直線」はどのようにして求めれば良いのでしょうか？

最小二乗法 (Ordinary Least Squares, OLS)

最小二乗法は、残差平方和RSSが最小になるように回帰係数（傾きと切片）を決定する手法です。

残差平方和 (Residual Sum of Squares, RSS)


図 1　残差

図 1にあるように残差を ε_{i}とすると、残差平方和RSSは次のようになります。

\def\arraystretch{1} \begin{align} RSS&=\sum_{i=1}^{n}(ε_i)^2 \tag{1}\\ &=\sum_{i=1}^{n}\{y_i-f(x_i)\}^2 \tag{2}\\ &=\sum_{i=1}^{n}\{y_i-(ax_i+b)\}^2 \tag{3} \end{align}

回帰係数

式(3)を a と b について偏微分したうえで式の変形を続けていくと、回帰係数が得られます。

\def\arraystretch{2} \begin{align} a&=\frac{\overline{xy}-\overline{x}\cdot\overline{y}}{\overline{x^2}-{\overline{x}}^2} \tag{4}\\ b&=\overline{y}-a\overline{x} \tag{5} \end{align}

式(4)については、実装のイメージをつかむために更に変形を進めていきます。それには、偏差平方和s_{xx}と偏差積和s_{xy}を使います。また、それぞれをデータ数nで割った分散σ_x^2と共分散σ_{xy}も使います。

\def\arraystretch{2} \begin{align} s_{xx}&=\sum_{i=1}^{n}(x_i-\overline{x})^2 \tag{6}\\ s_{xy}&=\sum_{i=1}^{n}(x_i-\overline{x})(y_i-\overline{y}) \tag{7}\\ σ_x^2&=\frac{1}{n}s_{xx}=\overline{x^2}-{\overline{x}}^2 \tag{8}\\ σ_{xy}&=\frac{1}{n}s_{xy}=\overline{xy}-\overline{x}\cdot\overline{y} \tag{9} \end{align}

式(4)の分子は共分散σ_{xy}そのもの、そして分母は分散σ_x^2そのものであることが分かります。

閉形式解

ここまでについて、偏差平方和s_{xx}や偏差積和s_{xy}との関係を合わせて整理して閉形式解としてまとめます。式(12)は式(5)の再掲です。

\def\arraystretch{2.5} \begin{align} a&=\frac{σ_{xy}}{σ_x^2}=\frac{s_{xy}}{s_{xx}} \tag{10}\\ &=\frac{\sum(x_i-\overline{x})(y_i-\overline{y})}{\sum(x_i-\overline{x})^2} \tag{11} \\ b&=\overline{y}-a\overline{x} \tag{12} \end{align}

実装

回帰係数を求めるols関数を作成して閉形式解を実装します。サンプルデータは、直線の方程式 y=1.6x+10 にノイズを加えて人工的に生成します。そのため面白みには欠けますが、この方程式に近い回帰係数が求まれば実装の正しさを確信できます。

Julia
using Statistics, CairoMakie, Random

# 最小二乗法（閉形式解）で回帰係数を計算
function ols(x::Vector, y::Vector)
    μx  = mean(x)                       # xの平均
    μy  = mean(y)                       # yの平均
    Sxy = sum((x .- μx) .* (y .- μy))   # 偏差積和
    Sxx = sum((x .- μx) .^ 2)           # 偏差平方和
    a   = Sxy / Sxx                     # 傾き
    b   = μy - a * μx                   # 切片
    (a = a, b = b)                      # 回帰係数
end

# サンプルデータを人工的に生成
function generate_sample_data()
    Random.seed!(0)
    x = collect(0:30)                   # xのデータを生成
    noise = randn(length(x)) * 4        # ノイズを生成（正規分布に従う）
    y = 1.6x .+ 10 .+ noise             # yのデータを生成（ノイズを付加）
    x, y
end

# サンプルデータと回帰直線を可視化
function vis_regression_line(x, y, coef)
    rl = coef.a .* x .+ coef.b          # 回帰直線
    a  = round(coef.a, digits=2)        # 傾き
    b  = round(coef.b, digits=2)        # 切片
    eq = "y = $(a)x + $(b)"             # 直線の方程式
    fig = Figure()
    ax = Axis(fig[1,1], xlabel="x", ylabel="y")
    scatter!(ax, x, y, color=:blue, markersize=10)
    lines!(ax, x, rl, color=:red, linewidth=2)
    text!(ax, 0, 56, text=eq, fontsize=20)
    ylims!(ax, 0, nothing)
    display(fig)
    fig
end

ハンズオン

Julia
x, y = generate_sample_data()           # サンプルデータを人工的に生成
coef = ols(x, y)                        # 回帰係数を取得
fig  = vis_regression_line(x, y, coef)  # サンプルデータと回帰直線を可視化
save("03-fig2.png", fig)


図 2　回帰直線

予測例

図 2の回帰直線 y=1.55x+10.54 を利用して、x=20 と x=40 について予測してみます。臨場感を持って予測を行うには「気温に対するビールの売上げを予測する」みたいなことを思い浮かべます。

\def\arraystretch{1.5} \begin{aligned} f(20)&=1.55・20+10.54=41.54\\ f(40)&=1.55・40+10.54=72.54 \end{aligned}

Discussion