エクセル近似曲線の罠 - 小人さんの妄想

エクセル近似曲線の「指数近似」「累乗近似」は、いわゆる非線形最小二乗法ではない。
　・エクセルで用いているのは、データの対数に直線をあてはめるという方法。
　・いわゆる非線形最小二乗法とは、残差（誤差）の２乗和を最小にする方法。

詳しいことは以下のページで尽きているのだが、エクセルの近似曲線は便利だと、
私自身多くの人に勧めている手前、注意を忘れないよう記載しておく。

* Excelグラフ累乗，指数，多項式近似の論文記載の注意(生物科学研究所井口研究室)
>> https://biolab.sakura.ne.jp/excel-graph.html

* 指数近似

赤い点で描かれているのが元になるデータ、
　Y = 0.8 * exp( 0.3 * X ) に、Xに比例する大きさで正規乱数を加えたもの。
下側の青い線がエクセルの指数近似曲線。
上側の緑の線はＲ言語による非線形回帰。

* 累乗近似

赤い点で描かれているのが元になるデータ、
　Y = 0.3 * X ^ 1.5 に正規乱数を加えたもの。
直線に近い青い線がエクセルの累乗近似曲線。
よりカーブのきつい、緑の線はＲ言語による非線形回帰。

ぱっと見に、あてはまりが良さそうなのはＲ言語の方に思える。
だからといってエクセルが“間違っている”というわけでもない。
エクセルはＲ言語と異なる方法で線を引いているのである。
Ｒ言語がもともとの回帰式の上で誤差を最小化するのに対し、
エクセルは両辺対数をとった式の上で誤差を最小化する。

* 指数近似
　　Y = a * X ^ b　　　-- もともとの回帰式
　　ln(Y) = ln(a) + b ln(X)　　　-- 両辺対数をとった式
* 累乗近似
　　Y = a * exp(b * X)　　　-- もともとの回帰式
　　ln(Y) = ln(a) + b * X　　　-- 両辺対数をとった式

もともとの回帰式は曲線だが、両辺対数をとった式は直線となる。
誤差の最小化は、直線の方がずっと簡単だ。
直線に直した上での回帰は以下のようになる。

* 指数近似

* 累乗近似

この直線で求めたパラメータが、エクセル近似曲線のものに一致していることが見て取れる。
　※ exp(0.2526)=1.2874, exp(-0.1613)=0.851
エクセル近似曲線は、直線回帰にちょっとオマケを付け足して実現した機能だったのだ。
エクセルの指数近似、累乗近似では、ゼロやマイナスを含むデータに線を引くことができない。
なぜできないかというと、ゼロやマイナスの対数がとれないからである。
また、エクセルには指数回帰を行う"LOGEST"という関数があるが、この関数で求めた数字も、近似曲線で求めた数字と同じである。
（もう１つの、直線回帰のための"LINEST"関数は特に問題無い。）

さらに、エクセルでは近似曲線に対して決定係数R2を計算することができるのだが、
決定係数R2は本来、線形回帰についてのあてはまり指標であり、非線形回帰に用いることはできない。
では、エクセルに出てくる決定係数は何なのかというと、対数をとって直線に直したときの決定係数だったのだ。
なので、この数字をうのみにして、そのまま用いるべきではない。
非線形回帰であてはまりの良さを評価するには、決定係数ではなく、
残差の２乗和、あるいは（回帰分析の）標準誤差を用いるのが妥当だ。

以上の懸念は、指数近似、累乗近似の話であって、他の種類の曲線「対数近似」「多項式近似」では心配しなくてよい。
実際試してみると、対数近似、多項式近似の結果はＲ言語と一致する。
なぜそうなるかというと、

　・対数近似，多項式近似は線形回帰
　・指数近似，累乗近似は非線形回帰
だからである。
形の上で、線形＝直線、非線形＝曲線、と思い込んでいる人もいるかもしれないが、
回帰分析の場合、回帰式が推定パラメータについて１次のものが線形、そうでないものが非線形である。
　・対数近似の回帰式：　y = a ln(X) + b
　・多項式近似の回帰式:　y = a X^2 + b X + c
いずれの回帰式も（Xではなく）a, b, c については１次式となっている。

で、指数近似，累乗近似がやりたかったとき、結局のところどうすれば良いのか。
・エクセルを使ったときは、エクセルの結果だと明記しておく。
・できればＲ言語、あるいは信頼できる統計用ソフトを使った方が良い。
・実はエクセルの「ソルバー」機能を使えば、非線形回帰も不可能ではない。
　でも、そこまでするなら素直にＲ言語使えば良いのではないかな。

* Excelのチャートオプション"近似曲線の追加"機能の評価
>> https://ci.nii.ac.jp/naid/110007025777

以下はＲ言語上での操作です。

### 累乗近似、非線形回帰
# データを用意
x <- c(1:12)
y <- c(
	1.075513124,
	2.104125184,
	2.045615008,
	3.70051273,
	2.221801362,
	2.769819828,
	6.656385201,
	7.037107115,
	7.976938733,
	8.150842204,
	11.70331758,
	13.66264136 )

nonline <- nls( y ~ a * x^b, ,start=c(a=1,b=1))  # 非線形回帰
summary(nonline)  # 結果の要約

Formula: y ~ a * x^b

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a   0.2680     0.1204   2.226   0.0502 .  
b   1.5587     0.1959   7.957 1.23e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 1.132 on 10 degrees of freedom

Number of iterations to convergence: 8 
Achieved convergence tolerance: 2.96e-06

# プロットしてみる
a <- coef(nonline)[1]
b <- coef(nonline)[2]
plot( y ~ x)
curve( a*x^b, add=T, col="red")

### 指数近似、非線形回帰
# データを用意
x <- c(1:12)
y <- c(
	3.283583359,
	2.754235941,
	1.674398192,
	3.556187361,
	0.976339451,
	4.652897765,
	15.21122786,
	12.09124224,
	8.109277666,
	13.59872153,
	23.046961,
	24.17609857 )

nonline <- nls( y ~ a * exp(b*x), ,start=c(a=1,b=1))  # 非線形回帰
summary(nonline)  # 結果の要約

Formula: y ~ a * exp(b * x)

Parameters:
  Estimate Std. Error t value Pr(>|t|)    
a  1.53717    0.67319   2.283 0.045521 *  
b  0.23269    0.04152   5.605 0.000226 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 3.404 on 10 degrees of freedom

Number of iterations to convergence: 15 
Achieved convergence tolerance: 3.005e-07

# プロットしてみる
a <- coef(nonline)[1]
b <- coef(nonline)[2]
plot( y ~ x)
curve(a * exp(b*x), add=T, col="red")