自動微分モジュールで遊ぶ
HackageDBでnumbersというパッケージを見つけました。この中にData.Number.Difという自動微分のモジュールが含まれていたのですが、これがなかなか面白いです。
サンプルコードを全然見かけないので、いろいろと書いてみます。
自動微分とは
自動微分というのは関数の導関数の値を自動的に求めるということですが、これには「数値的」にやる方法と「記号的」にやる方法があります。数値的にというのは
という式のhに適当な微小量を与えて、近似的に導関数の値を求める方法を言います。
一方記号的にというのは
の様に、微分公式を使用して具体的に関数を求める方法を言います。
Difモジュールが行うのは後者ですが、Haskellの遅延評価と多相型を生かして非常に使い勝手の良いものになっています。プリプロセッサを通すとかいった手間もありません。
自動微分のアルゴリズムについてはAutomatic differentiationを参照して下さい。
簡単な例
DifモジュールではHaskellの数値型の多相性を利用して、通常どおりに定義した関数をそのまま微分できるようにしています。
通常どおりというのは、数式に限らずif文を利用した関数なども含みます。
> :m +Data.Number.Dif > let f(x) = x^2 + 3*x + 2 > f 2 12.0 > deriv f 2 7.0 > deriv (deriv f) 2 2.0
上ではの微分をderivという関数によって行っています。
最初のf 2の部分は通常の関数と通常の数値ですが、deriv f 2の引数のfは微分を計算するために、Dif a -> Dif aという特殊な型になっています。
ここで、Dif aという型がNumクラスやFloating等のすべての数値クラスに属しているので、通常の式と全く同じく記述することができます。Haskellでは数値リテラルも多相的に振る舞ってくれます。
また、sinなどの無限回微分出来る関数も遅延評価を利用して自然に表現されています。
応用編1
ニュートン法です。
import Data.Number.Dif newton f init eps = let f' = deriv f next = init - (unDif f init) / (f' init) in if abs (next - init) < eps then init else newton f next eps
関数fと初期値、許容誤差を与えるとf(x) = 0の解の近似値が求まります。
*2
> let e = 1.0e-15 > newton (\x -> x**2 - 2) 1 e 1.4142135623730951 > newton (\x -> x**2 - 3) 1 e 1.7320508075688772 > newton sin 3 e 3.141592653589793 > newton (\x -> log x - 1) 2 e 2.7182818284590455
応用編2
Gnuplotと連携させてグラフを書いてみます。例として導関数と接線を書いてみます。f(x) = ...の部分を書き換えればいろいろな関数の微分が試せます。
import Data.Number.Dif import Data.List import System.IO import System.Cmd -- x=aのまわりでの一次近似 (接線) tangent f a = \x -> (deriv f a) * (x - a) + (unDif f a) main = let -- 変域設定 xmin = -2*pi xmax = 2*pi step = 0.1 -- 表示する関数の用意 f(x) = sin x funs = [f, deriv f, tangent f 1.0, tangent f 2.0] in let range = [xmin, xmin+step .. xmax] graphs = [map f range| f <- funs] datas = transpose (range : graphs) in do -- データの用意 (datafile, hd) <- openTempFile "/tmp" "diff-data" mapM (hPutStrLn hd) $ [concat . intersperse "\t" . map show $ d | d <- datas] hFlush hd -- バッチファイルの用意 (cmdfile, hc) <- openTempFile "/tmp" "diff-command" hPutStrLn hc"set xzeroaxis" hPutStrLn hc"set yzeroaxis" hPutStrLn hc $ "plot " ++ ( concat . intersperse "," . map (\n -> show datafile ++ " using 1:"++show n++" w l") ) [2..length graphs + 1] hPutStrLn hc "pause -1" hFlush hc system $ "gnuplot " ++ cmdfile
