重み付き最小二乗法を用いた対数関数の回帰分析

はじめに
通常の最小二乗法ではだめか？
重み付き最小二乗法
Rを用いた重み付き最小二乗法
Box-Cox変換
まとめ
参考書籍

はじめに

自分が大学生だったころ、

$y=e^{a+bx}$

の形（C=e^aとすればy=Ce^bxの形）で最小二乗法による回帰分析を行う機会が結構ありました。

当時は、両辺logをとって

$Y=\log y=a+bx$

と変換してから最小二乗法を適用すればいいと教わったのですが、本当にそれでいいのかな～と少し気になっていたので、本記事を書くことにしました。

通常の最小二乗法ではだめか？

y=e^a+bxの両辺logをとって単純に最小二乗法を適用することに、なんとなく気持ち悪さを感じて調べてみたら、案の定分散σ²に関して注意する必要があることが分かりました。

通常の最小二乗法では、すべての実験値yが等しい独立な分散σ²の正規分布に従うことを仮定しています。そのため、対数変換すると、分散σ²も対数変換されるので等分散の仮定が破綻してしまうんですね。

母数yの区間推定値を実験値y_iと分散σ²（標準偏差σ）を用いて

$y=y_i\pm\sigma$

と表すと、logをとれば

$\displaystyle{ Y=log(y)=log(y_i\pm\sigma) \approx log(y_i)\pm\frac{\sigma}{y_i}}$

であることが分かります。つまり、もとのσが一定だとすると、Y_iの分散はσ/y_iとなり実験値y_iの大きさによって異なってしまうことがわかります。

したがって、変化後のY=a+bxに対して、σ一定の下でいつもの最小二乗法を適用することは原理的には間違っているということになるでしょう。

※ただし、yの分散が相対的、つまりσがyに比例して大きくなる場合にはσ_Y=σ/yが一定となるため、logで変換した後に通常の最小二乗法を適用できます。

重み付き最小二乗法

σがyの値によって異なる場合には、重み付き最小二乗法を使うのが推奨されます。

重み付き最小二乗法の式は、最尤法に基づいて最小二乗法を導出するとき自然に現れるため、実は重み付きの方が通常の最小二乗法より一般的な方法になります。

xが母平均µの正規分布に従って観測されるとき、y=y₁, y₂を観測する確率は

$\displaystyle{ P(y_1, y_2)=\frac{1}{2\pi\sigma_1\sigma_2}exp\left(-\frac{(y_1-a-bx_1)^2}{\sigma_1^2}-\frac{(y_2-a-bx_2)^2}{\sigma_2^2}\right) }$

となります。ここで、σ₁、σ₂はそれぞれy₁、y₂の標準偏差です。最尤法の考えではこの確率を最大化するので、

$\displaystyle {\frac{(y_1-a-bx_1)^2}{\sigma_1^2}+\frac{(y_2-a-bx_2)^2}{\sigma_2^2}}$

を最小化するaとbを求めればよいことになります。yが2つだけでなく、n個ある場合には

$\displaystyle{ min \sum^n_i{\frac{(y_i-a-bx_i)^2}{\sigma_i^2}} }$

となります（σ一定とすると最小二乗法と同じ意味になりますね）。

この式より、σを一定としない場合には通常の最小二乗法に重み1/σ_i²をかければよいことが分かりました。したがって、

■ポイント■

logy=Y=a+bxの重み付き最小二乗法における重みは

$\displaystyle{ weight=\frac{1}{\sigma_{Yi}^2}=\frac{1}{(\sigma_i/y_i)^2}=\frac{y_i^2}{\sigma_i^2} }$

と設定すればよいことになります。

Rを用いた重み付き最小二乗法

Rを用いた重み付き最小二乗法の例を示します。以下コードではyの分散をs2=1E-15で一定としました。logyはx=0におけるyの値y0で正規化した（1になるようにした）値です。

そして、重み付き最小二乗法と重みなしの通常の最小二乗法によってlogy~xの回帰分析を行っています。

x <- c(0, 1, 2, 3, 4, 5)
y <- c(5.7E-06,2.8E-06,1.1E-06,5E-07,3E-07,7.0E-08)
s2 <- 1E-15
#log変換してプロット
logy <- log(y/y[1])
plot(x, logy, xlab = "x", ylab = "log(y/y0)", pch=19, ylim=c(-5,0))
#重みの計算
w <- y/sqrt(s2)
#重み付き最小二乗法
result <- lm(logy~x, weights =w^2)
#通常の最小二乗法
result2 <- lm(logy~x)
abline(result, col="blue")
abline(result2, col="red")
#エラーバー。Yの分散(yに対してはs2で一定だが、Yに対してはy/sqrt(s2)
arrows(x, logy - 1/w, x, logy + 1/w, angle=90, code=3, length = 0.1)

図1に、logy vs. xのプロットを示します。yに対する分散s2は一定ですが、logyのプロット中のエラーバーはs2/yなので、観測点によって大きさが異なっています。

また、青線が重み付き最小二乗法による回帰直線、赤線が通常の最小二乗法による回帰直線です。重み付き最小二乗法（青）の方が、yの大きな点（エラーバーの小さな点）によりフィットした線になっています。一方、通常の最小二乗法（赤）はエラーバーの幅に関係なく、すべての点にフィットするような直線ですね。

図1 logy vs. xのプロット。重み付き最小二乗法(青)は通常の最小二乗法(赤)と異なり、誤差の大きなx=5の点の影響をあまり受けていないことが分かる。

ちなみに、図1では重み付きと重みなしの最小二乗法で明らかに異なる結果になりましたが、実験点がきれいに一直線上にある場合には、どちらの最小二乗法で回帰分析を行っても実用上の問題はないようです。

Box-Cox変換

回帰分析で目的変数の誤差が正規分布に従わない場合、目的変数を適当に変換して、最小二乗法を正しく実行できるようにする、という今回の内容は、実はBox-Cox変換という話題につながっています。もしご興味ありましたら、こちらの記事もご覧ください。

sturgeon.hatenablog.com

まとめ

本記事では、対数に関する回帰分析の際に用いる重み付き最小二乗法についてまとめました。y=e^a+bxの形で表されるyに対して分散σ²が一定の場合、logyと変換して通常の最小二乗法を適用すると、変換後は分散が一定ではなく(σ/y)²になるため、重み付き最小二乗法を使います。重みは(y/σ)²と設定します。