乱数降下法(1) 関数推定
いくつかの変数の組があって、それらの間の関係式の形がわかっているとき、式の係数を探したいときがある。たとえば、次のようなデータがあるとする。
#x y z f(x,y)=z 10 3 -19 16 9 -19 9 2 -19 15 16 -19 12 12 -19 17 16 -19 5 13 -17 4 19 -15 15 3 -19 13 19 -19 16 14 -19 3 10 -10 17 10 -19 18 20 -19 12 16 -19 4 15 -15 18 13 -19 1 20 16 8 9 -19 8 7 -19
このとき、関数f(x,y)が単純な関数だったり、微分可能な関数であれば手計算で求められるが、面倒な形をしていることも多い。例の関数は次の形をしている。
(defun f (a b p) (lambda (x y) (floor (- (* a (expt p x)) b))))
このa,b,pに乱数を与えて例題を作ってみた。係数は、それぞれ範囲は分かっているが、実際の数値は不明なものとしている。
(let nil (setq arg (list (1+ (random 100)) (1+ (random 100)) (random 2.0))) nil) (setq test-case (mapcar (lambda (x) (append x (list (apply (apply 'f arg) x)))) (mapcar (lambda (_) (list (1+ (random 20)) (1+ (random 20)))) (make-list 20))))
test-caseは、上の内容が*1。
(defun test (a b p) (lambda (x y z) (- z (floor (- (* a (expt p x)) b)))))
a,b,pが正しければ解が0になるはずである。解の範囲はa,b,q->1-100,1-100,0-2なので、乱数降下法の関数rdmに関数、テストケースとともに与える。
(rdm 'test test-case '((1 100)(1 100) (0 2))) |(9 (84.20693 18.02192 0.4090839)) |(4 (63.83097 18.40089 0.5271799)) |(10 (92.37186 18.87796 0.4470151))
実行すると、導出した係数と、それをtest-caseに当てはめたときの最小値との距離が返ってくる。指数を含む3変数のため、最小値(距離=0)はすぐには見つけられない。何度か実行して変数の範囲を見て、指定し直す。
(rdm 'test test-case '((40 100)(1 50) (0 1))) |(1 (72.83604 18.4344 0.4788978)) |(3 (67.85478 18.36407 0.4844662)) |(1 (73.70525 18.48913 0.4906388))
すると、また距離と範囲が小さくなってるのが分かるので、範囲をより狭くする。そろそろ解が見つかりそうなので、最小値以外の結果を弾くように指定する。
(rdm 'test test-case '((60 80)(10 27) (0.3 0.6)) :valid 0) |(0 (69.88915 18.74341 0.5080137)) |(0 (69.82953 18.94899 0.5098826)) |(0 (68.81656 18.76009 0.507199))
解がいくつか出てくるが、問題の関数はfloorによって切り下げをしているため、適合する係数に幅がある。そこで何点かサンプルをとって範囲を調べる。
(rdm 'test test-case '((60 80)(10 27) (0.3 0.6)) :valid 0 :sample 30) |((67.46078 72.46016) (18.23532 18.96503) (0.4854829 0.5193323))
結果はこういう感じ。切り下げだから、整数部分はもう少し上に幅広く見た方が良いかも知れない。まだ範囲が広めだが、結果がステップ状になっている関数に対してテストケースが少ないので仕方がない。答え合わせをすると、
arg |(70 19 0.5123419)
なので、ちゃんとそれらしい結果は出ている。
この乱数降下法の関数は以下のものである。結構処理が重いので(compile 'rdm)など、コンパイルしてから使った方がよい。
(defun rdm (test test-case range &key mode verbose (rate 1) valid sample) (and valid (<= valid 0) (setq valid 1e-6)) (when (or (null mode) (not (listp mode))) (setq mode (let* ((f (if (numberp mode) (max 1 mode) (length range))) (l (append (make-list (min 8 (1- f)) :initial-element 2) (list (ceiling (* 22 (expt 1.35 f))))))) (push (max 2 (- 26 (* 2 f))) l)))) (when verbose (format t "mode:~A/acc:~A/amount:~A~%" mode (* rate (expt 0.5 (- (apply '+ (cdr (reverse mode))) (length mode)))) (1+ (apply '* mode))) (refresh-screen)) (let ((show (1- (length mode))) (time (get-internal-real-time)) (gen (lambda (x) (+ (car x) (random (float (- (cadr x) (car x)))))))) (labels ((sample (range) (let ((args (mapcar gen range))) (list (apply '+ (mapcar (lambda (x) (abs (apply (apply test args) x))) test-case)) args))) (narrow (range s conv) (mapcar (lambda (x y) (let ((lim (* conv (/ (- (cadr x) (car x)) 2)))) (list (max (car x) (- y lim)) (min (cadr x) (+ y lim))))) range (cadr s))) (app-sub (model range conv limit) (let ((rest (car limit))) (if (cdr limit) (progn (while (plusp (1+ (decf rest))) (let ((new (app-sub model range conv (cdr limit)))) (and new (> (car model) (car new)) (setq model new)) (setq conv (* conv 0.5)))) (when (and verbose (= (length limit) show)) (format t ">~A~%" model) (refresh-screen)) (keep-response 500)) (while (plusp (1+ (decf rest))) (let ((new (sample (narrow range model conv)))) (and (> (car model) (car new)) (setq model new)))))) model) (keep-response (interval) (when (< (+ time interval) (get-internal-real-time)) (do-events) (setq time (get-internal-real-time))))) (long-operation (if (numberp valid) (let ((r (mapcar (lambda (x) (list (cadr x) (car x))) range))) (loop (let ((res #1=(app-sub (sample range) range rate mode))) (when (<= (abs (car res)) (abs valid)) (if (null sample) (return res)) (setq r (mapcar (lambda (x y) (list (min y (car x)) (max y (cadr x)))) r (cadr res))) (if (< (decf sample) 0) (return r)))))) #1#)))))
*1:x y z)...)という形で格納されている。今回はこれを元に、乱数降下法を用いて係数を探してみる。 乱数降下法は、解の存在が期待される空間内で乱数を取り、解に近い位置へ探索範囲を狭めて解を探す手法である。ちゃんと調べたわけじゃないので名称が間違ってるかも知れないが、動作には関係ない。基本的には降下法による探索なので、解の周辺が解の方向へ勾配していることが前提だが、多少起伏があったり、微分不可能な区間があっても探索できるのが強みである。欠点は、収束が遅いことと、収束のクセによって局所最小値に嵌ったり範囲の端付近にある解を見つけにくいことである。 とりあえず実際に試してみよう。まず、関数を乱数降下法に適した形に直す。正しい係数が関数に与えられたときに最小の値を取ればよいから、次のような形にする((乱数降下法の関数がf(x,y)=zの形で評価しても良いが、f(x,y)-z=0の方が一般性がある。
