はじめよう実験計画

実験を早く終わらせるための技術

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

はじめに

通常の回帰分析はy=ax+bの形で行いますが、以前

 y=e^{a+bx}

の形(C=eaとすればy=Cebxの形)で最小二乗法による回帰分析を行う機会がありました。

その時、両辺logをとって Y=logy=a+bx に変換してから最小二乗法を適用すればいいと教わったのですが、本当にそれでいいのか少し気になっていたので、本記事を書くことにしました。

 

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

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

通常の最小二乗法では、すべての実験値yが等しい独立な分散σ2正規分布に従うことを仮定しています。そのため、y=ea+bxをY(=logy)=a+bxに変換すると、分散σ2も変換されて等分散の仮定が破綻します。

母数yの区間推定値を実験値yiと分散σ2標準偏差σ)を用いて

 y=y_i\pm\sigma

と表すと、logをとれば

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

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

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

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

 

重み付き最小二乗法

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

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

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

 \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) }

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

 \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/σi2をかければよいことが分かりました。

したがって、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=ea+bxの形で表されるyに対して分散σ2が一定の場合、logyと変換して通常の最小二乗法を適用すると、変換後は分散が一定ではなく(σ/y)2になるため、重み付き最小二乗法を使います。重みは(y/σ)2と設定します。

参考書籍

最小二乗法や誤差に関して、こちらの本が非常に役に立つと思います。この本には大学時代、実験課題のたびにお世話になりました。誤差に関してどのように議論したらいいのか、結構それだけでも考察になるんですよね。