自分のコピペ用にPerlで誤差関数と逆誤差関数を書いた

学習器を実装しようとすると唐突に逆誤差関数(erf(誤差関数)の逆関数)が必要になったりする。
こんなときに慌てず騒がずコピペできるようにPerlで実装したものをメモしておく。

誤差関数および逆誤差関数は解析的に解が求まらないけれどテイラー展開で近似することができる(参考: http://ja.wikipedia.org/wiki/%E8%AA%A4%E5%B7%AE%E9%96%A2%E6%95%B0)。

#!/usr/bin/perl
use strict;
use warnings;
use Math::Trig qw/pi/;

my $num = shift @ARGV;
my $pre  = 0;
for my $n (1..10) {
    my $order   = $n * 10;
    my $erf_inv = erf_inv($num, $order);
    my $diff    = abs($erf_inv - $pre);
    print "$order : $erf_inv($diff)\n";
    $pre = $erf_inv;
}


sub erf {
    my ($z, $order) = @_;

    my $value = 1;
    my $term  = 1;
    for my $n (1 .. $order) {
        $term  *= (-1 * $z * $z / $n);
        $value += ($term / (2 * $n + 1));
    }
    return (2 * $z * $value / sqrt(pi()));
}

sub erf_inv {
    my ($z, $order) = @_;

    my $value  = 1;
    my $term   = 1;
    my @c_memo = (1);
    for my $n (1 .. $order) {
        $term  *= (pi() * $z * $z / 4);
        my $c = 0;
        for my $m (0 .. ($n - 1)) {
            $c += ($c_memo[$m] * $c_memo[$n - 1 - $m] /
                   ($m + 1) / (2 * $m + 1));
        }
        push(@c_memo, $c);
        $value += ($c * $term / (2 * $n + 1));
    }
    return (sqrt(pi()) * $z * $value / 2);
}

これをtest.plという名前で保存して実行する。絶対値が1に近い値だと次数を大きくしても誤差が残る。

$$ perl test.pl 0.5
10 : 0.476936272711668(0.476936272711668)
20 : 0.476936276204468(3.49279971612049e-09)
30 : 0.47693627620447(1.66533453693773e-15)
40 : 0.47693627620447(0)
50 : 0.47693627620447(0)
60 : 0.47693627620447(0)
70 : 0.47693627620447(0)
80 : 0.47693627620447(0)
90 : 0.47693627620447(0)
100 : 0.47693627620447(0)

$$ perl test.pl 0.9
10 : 1.15506001297562(1.15506001297562)
20 : 1.16254511208024(0.00748509910461626)
30 : 1.16304179407171(0.000496681991466019)
40 : 1.16308296621338(4.11721416719857e-05)
50 : 1.1630867444565(3.77824312502995e-06)
60 : 1.1630871121658(3.67709294790686e-07)
70 : 1.16308714935333(3.71875288340817e-08)
80 : 1.16308715321746(3.86413634423377e-09)
90 : 1.16308715362715(4.09687617164423e-10)
100 : 1.16308715367127(4.41173764187397e-11)