廿TT

譬如水怙牛過窓櫺 頭角四蹄都過了 因甚麼尾巴過不得

一次元の解探索、f(x)=0を解く、二分法(bisection method)

大雑把にいうと、


f(a)とf(b)の符号が違えば、aとbの間にf(c)=0となるようなcがあるはずだ
だから区間[a,b]を半分づつ狭めていってcに近づこう


というのが二分法の考え方である。
二分法(wikipedia)


なんてわかりやすいアルゴリズムなんだ
こういうのがいい。複雑なことはわからないんだよ
二分法が精神的支柱


以下で、
Maria L. Rizzo "Statistical Computing with R"の
Example 11.6(f(x)=0を解く)を見ていく

Statistical Computing with R (Chapman & Hall/CRC The R Series)

Statistical Computing with R (Chapman & Hall/CRC The R Series)

 f(y)=a^2 + y^2 + \frac{2ay}{n-1}-(n-2)=0
ここでaは指定された定数、nは2以上の整数
これをyについて解く
今回は
a=1/2,n=20
とする
以下で正の解を探す

f <- function(y, a, n) {
		a^2 + y^2 + 2*a*y/(n-1) - (n-2)
}  #まず関数を定義した
a <- 0.5
n <- 20
b0 <- 0  
b1 <- 5*n  #最初の区間を[0,5n]とした

#二分法を使う
it <- 0  #繰り返しの回数を数えるのに使うことにする
eps <- .Machine$double.eps^0.25  # 許容される誤差の範囲を ".Machine$double.eps^0.25" の値にする
r <- seq(b0, b1, length=3)		# r=c(0,50,100)
y <- c(f(r[1], a, n), f(r[2], a, n), f(r[3], a, n))
if (y[1] * y[3] > 0)stop("f does not have opposite sign at endpoints")
#y[1]とy[3]が同符号ならばstop
#stopは、現在の式の実行を中断し、引数として与えられたメッセージを表示する関数

while(it < 1000 && abs(y[2]) > eps){
		it <- it + 1
		if (y[1]*y[2] < 0) {	# y[1]とy[2]が異符号ならば,
			r[3] <- r[2]	# r[2]を区間の右端に更新
			y[3] <- y[2]
		} else {			#そうでなければ,
			r[1] <- r[2]	# r[2]を区間の左端に更新
			y[1] <- y[2]
		}
		r[2] <- (r[1] + r[3])/2  #r[2]は区間の中点
		y[2] <- f(r[2], a=a, n=n)
		print(c(r[1], y[1], y[3]-y[2]))
}

これを実行すると

[1]    0.000  -17.750 1876.316
[1]   0.0000 -17.7500 469.4079
[1]   0.0000 -17.7500 117.5164
[1]   0.00000 -17.75000  29.46135
[1]  3.125000 -7.819901 17.172081
[1]  3.125000 -7.819901  6.754986
[1]  3.906250 -2.285619  3.530081
[1]  3.906250 -2.285619  1.650599
[1]  4.1015625 -0.7113133  0.8348365
[1]  4.1015625 -0.7113133  0.4102657
[1]  4.1503906 -0.3058160  0.2057289
[1]  4.1748047 -0.1012793  0.1030135
[1]  4.17480469 -0.10127926  0.05139497
[1]  4.18090820 -0.04995880  0.02570680
[1]  4.18395996 -0.02427063  0.01285573
[1]  4.185485840 -0.011419556  0.006428445
[1]  4.186248779 -0.004992275  0.003214368
[1]  4.186630249 -0.001778197  0.001607221
[1]  4.1868209839 -0.0001710497  0.0008036194
[1]  4.1868209839 -0.0001710497  0.0004018029
[1]  4.1868209839 -0.0001710497  0.0002008997

なんかいっぱいでた
結果は

> r[2]  # r[2]は区間の中点
[1] 4.186845
> y[2]  # y[2]はf(r[2])の値
[1] 2.984885e-05
> it  #繰り返しの回数
[1] 21

21回の繰り返しで
x=4.186845のときf(x)=0.00002984885
であることがわかった。

許容される誤差の範囲にした.Machine$double.eps^0.25は、ぼくの機械だと0.0001220703なので
0.00002984885 < 0.0001220703
になっていることがわかる


しかし、これがどのくらい良い近似なのかとかはよくわからないね
.Machine$double.epsの意味もわからない


この方程式はふつうに(二次方程式の解の公式で)解くと、
 y=\frac{-a}{n-1} \pm \sqrt{n-2 - a^2+(\frac{a}{n-1})^2}


テキストだと
 y=\frac{-a}{n-1} \pm \sqrt{n-2 + a^2+(\frac{a}{n-1})^2}
になっているが、ただの間違いであろう


なので正確な解は

> c(-a/(n-1)+sqrt(n-2 - a^2 + (a/(n-1))^2),-a/(n-1)-sqrt(n-2 - a^2 + (a/(n-1))^2))
[1] 4.186841 -4.239473

である。


次はブレント法なんだけど、どうなの?大丈夫?いけるのか?おれ
ブレント法(wikipedia)