乱数降下法(2) 関数を解く
今回も乱数降下法(1) 関数推定 - 象徴ヶ淵で作ったrdmの話である。係数決定に使えるのだから、当然普通に方程式を解くことも出来る。
(defun f1 (x) (lambda () (+ 11 (* -3 x) (* x x))))
f1のように、変数を持たない関数を返すようにすれば、その式の結果の絶対値が最小、つまりf1(x)=0になるようなxの値を求めてくれる。この場合のテストケースは空のリストを1つ入れておけばよい。
(rdm 'f1 '(()) '((-7 7))) |(8.75 (1.500338))
実際に試してみると、絶対値が0にならず、解を求めることが出来ない。f1ではxは実数しか取れないので解くことが出来ないのだ。こういうときは複素数空間で探すように書き換える必要がある。
(defun f2 (a b) (lambda () (let ((x (complex a b))) (+ 11 (* -3 x) (* x x))))) (rdm 'f2 '(()) '((-7 7) (-7 7))) |(1.907349e-6 (1.5 -2.95804)) |(2.132481e-6 (1.5 2.95804)) (funcall (f2 1.5 -2.95804)) |#C(-9.536743e-7 0.0) (funcall (f2 1.5 2.95804)) |#C(-9.536743e-7 0.0)
距離が0になっていないが、数字にe-6が付いてるのは桁が6桁小さいことを表しているので、ほぼ0とみなして良い。これは単に処理系の精度の問題である。解が2つあるので、何度か試すと2種類出てくるはずだ。一度の探索で1つの解しか見つけないが、適当な確率で解のいずれかに収束するので、何度も繰り返せば全ての解を探すことが出来る。
(defun f3 (a) (lambda () (* (- a 5) (- a 4) (- a 3) (- a 2) (- a 1) a (- a -1)(- a -2)(- a -3) (- a -4) (- a -5))))
上のように、解が多数有る場合を考える。複雑ではあるが、結局1変数関数なのでrdmにとっては問題にならない。
(rdm 'f3 '(()) '((-10 10))) |(0.0 (-3.0)) |(0.0 (-1.0)) |(0.0 (2.0)) |(0.0 (-1.0))
ほとんどの試行でスコアが0に収束する。ここでは解の範囲を-10〜10としたが、範囲を広げても大体問題ない*1。さて、数回実行した程度では全ての解が見つからないので、繰り返し探索して解が全て求まるか試してみよう。
その前に、デフォルトの設定では計算に時間がかかるので調整をしよう。キーのverboseをtにすると詳細情報が表示される。
(rdm 'f3 '(()) '((-10 10)) :verbose t) |mode:(24 30)/acc:2.384186e-7/amount:721 |(0.0 (-1.0))
ここでmodeというのは探索方法で、要素の個数が階数、数値が繰り返し数を表す。accは何桁まで計算するか、amountは試行回数を示す。accは大雑把にmodeの最初の数値に依存して、modeの数値が1減るとaccが倍になる。試行回数はほぼmodeの数値を掛け合わせた数である。今は出来るだけ厳密に解を求めたいのではなく、探索結果の分布を調べたいのでできるだけ少ない試行回数にしたい。
今回はいずれかの整数の近傍に収まりさえすればいいので、せいぜい3桁程度の精度で十分だ。modeの変更はキーのmodeで行う。
(rdm 'f3 '(()) '((-10 10)) :verbose t :mode '(20 5)) |mode:(20 5)/acc:3.814697e-6/amount:101 |(0.03965374 (-1.999999)) (rdm 'f3 '(()) '((-10 10)) :verbose t :mode '(16 3)) |mode:(16 3)/acc:6.103516e-5/amount:49 |(0.6571106 (0.999962))
上の結果だと50回か、もっと少なくても良いように見えるが、試行回数が少ないほど最小値に収束しない確率が増えるので、少し多めにとって100回くらいのものにする。
小数値だとブレがあって数えにくいので、解に近いかどうかで数える。解から離れているものは弾く。試行回数を数えるスクリプトは次のような感じ。割と遅いので試行回数は100k回程度とした。
(setq x 0 y 0 lst (make-list 11 :initial-element 0)) (loop (case x (1000 #1=(format t "amount ~A fail ~A~%~{~@6A~}~%~{~@6A~}~%~%" x y (list -5 -4 -3 -2 -1 0 1 2 3 4 5) lst) (refresh-screen)) (10000 #1# (refresh-screen)) (40000 #1# (refresh-screen)) (100000 #1# (return))) (do-events) (incf x) (multiple-value-bind (r a) (round (caadr (rdm 'f3 '(()) '((-10 10)) :mode '(20 5)))) (if (> (abs a) 0.001) (incf y) (if (<= 0 (incf r 5) 10) (incf (elt lst r)) (incf y)))))
結果。
amount 1000 fail 1
-5 -4 -3 -2 -1 0 1 2 3 4 5
0 15 48 94 219 262 198 121 35 7 0
amount 10000 fail 5
-5 -4 -3 -2 -1 0 1 2 3 4 5
4 85 401 1102 2040 2620 2077 1185 403 73 5
amount 40000 fail 26
-5 -4 -3 -2 -1 0 1 2 3 4 5
33 335 1568 4421 8334 10470 8287 4525 1639 332 30
amount 100000 fail 51
-5 -4 -3 -2 -1 0 1 2 3 4 5
94 798 3898 11064 20903 26105 20886 11329 3959 836 77±5になる結果が0になる結果の0.1％もなく、失敗に埋もれてしまいそうなくらいである。偏り方は正規分布に近く、端の方ほど見つかりにくい。探索範囲に入っている期間が長いほど補足される可能性が高いので、近くに他の解が多い解に流れる可能性が高くなる。サイコロを2個振ったとき7が一番多く出るのと同じようなものだ。この現象は乱数を複数回使っている限り逃れることは出来ない。
