Gaussian Process で滑らかな関数を発生させる

Gaussian Process は関数空間上の確率変数の一種。

#! /usr/bin/env perl
use strict;
use warnings;
use Math::Random qw/:all/;

my @xs = map { ($_/30) } (0..20);  # 関数への入力値列
my $n = scalar @xs;
my @v = map { my $x=$_; [map {exp -1/2*( ($x-$_)**2 )} @xs]} @xs; # カーネル
#my @v = map { [map 1, @xs] } @xs;
my @m = ((0)x$n);
print map join(qq{\n}, @$_).qq{\n\n},
  random_multivariate_normal(2, @m, @v); # 2個の出力値列が得られる。各列は入力値列に対応する。

入力値の刻みを細かくしたり、入力値の数を増やそうとすると破綻する。
たぶん、浮動小数点の誤差が問題になっている。

もっと細かく値をとりたいときは、
逆関数法で random_multivariate_normal
を実装するか、R使っとけと。

use bignum して exp を 0 周辺でテイラー展開するとか、
covariance matrix が対称行列になるように対象位置の値はコピーで作るとか、
やってみたが無効だった。

GPML pp.13--14 辺りを見て実装した。
# ちなみに入力を多次元化するときは、カーネルの$x-$_に square root 距離をかませればいいらしい。
GP が滑らかな関数を発生させる理由は、
カーネル関数の値、つまり入力値同士の「類似度」を covariance に指定しているから。
遠い入力値に対応する出力値同士は covariance が小さいので相関がないが、
近い場合は相関が強くて、
結果として、滑らかっぽい関数が出来上がる。
「前後の点の座標値とあまり離れない」という意味での滑らかさではなく、
「入力値が近ければ出力値が近い」という意味