PerlとJAGSで赤池情報量基準(AIC): AICの計算

[17/08/04 18:37 更新]

AIC等の計算をするサンプルスクリプトはtest.plというファイル名で用意してある。このスクリプトは、cubic00.csvのデータをpoly_model.jagのモデルで推定する。

データ(cubic00.csv)

ここで使うデータは、mkdata_00.plで生成した。3次関数に正規乱数を足したものだ。実行にはPDLが必要である。パラメータを変えて実行すれば、様子が違うデータが得られるはずだ。ここでは、alpha = [0, 1, 0, 0.1]、sigma = 2(tau=0.25)とした。ちょっと見には、直線なのか高次の項が入ってるのか迷うようなデータになっている。

cubic.png

モデル(poly_model.jag)

今回のモデルは、こんな感じ。 


model{
	x.bar <- mean(x)
	for (i in 1:N) {
		for (j in 1:K) {
			xij[i, j]  <- (x[i] - x.bar)^(j-1)
		}
	}
	for (i in 1:N) {
		y[i] ~ dnorm(mu[i], tau)
	}
	mu <- xij %*% alpha

	for (i in 1:K) {
		  alpha[i] ~ dnorm(0, 1.0E-6)
    }

	tau ~ dgamma(1.0E-6, 1.0E-6)
	sigma2 <- 1/tau
}

ちょっと分かりにくいかもしれないが、(x[i] - x.bar)の多項式の周りに正規分布しているというモデルだ。ここで、x.barはx[i]の平均値。おそらく、こうした方が、x[i]の多項式にするより収束がよい。

サンプルスクリプト(test.pl)

以下に、サンプルスクリプトtest.plをもとに、本パッケージの使い方を説明する。

スクリプトの前半は、多項式の項数Kを2から6まで振って推定してること以外は、reg_test.plとあまり変わらない。[ちなみに、項数がKだと、次数はK-1である] $coda->summaryのあとで、 


	$param = $jk->param;
	@alpha = @{$param->get("alpha")->value};
	push @alpha, 0;
	$tau = $param->get("tau")->value;

としているのは、今回の推定値を次回の初期値に設定するための工夫である。

そのあとで、AIC等を計算し、結果を配列に収めている。 


	my $aic = $jk->aic;
	my $aicc = $jk->aicc;
	my $bic = $jk->bic;
	push @result, "K = $k, AIC = $aic, AICc = $aicc, BIC = $bic\n";

ループを抜けてから、結果の表示。 


print @result;

ここでmainは終わりなのだが、次にMakeshift::Jagsを拡張している。これは、モデルに依存している部分なので、Makeshift/Jags.pmには書きたくない。 


package Makeshift::Jags;

use PDL ();
use PDL::Constants qw(PI);
use PDL::GSLSF::POLY;

# ic_param: give parameters to calculate xIC
# return: hash ref {loglik => log likelihood, nparam => number of parameters,
#                   ndata => number of data points}
sub xic_param {
	my $self = shift;
	my $x = PDL->new($self->data->get("x")->{value});
	my $xbar = $x->avg;
	my $y = PDL->new($self->data->get("y")->{value});
	my $alpha = PDL->new($self->param->get("alpha")->value);
	my $tau = $self->param->get("tau")->value;
	my $mu = gsl_poly_eval($x - $xbar, $alpha);
	my $lik = sqrt($tau / (2*PI)) * exp(-$tau*($y - $mu)**2/2);
	my $loglik = log($lik)->sum;
	my $ret = {loglik => $loglik, nparam => $alpha->nelem + 1,
			   ndata => $x->nelem};

	return $ret;
}

ここではPDLを使っている。そうでないと、もっと長いコードになるうえに、計算速度も落ちるだろう。

これは、対数尤度、パラメータ数、データ点数をハッシュに収めて返すルーチンである。サブルーチン名(xic_param)と、ハッシュのキー(loglik、nparam、ndata)は変えないでほしい。あとの部分は、どんな名前でどうやって計算してもよい。ちなみに、パラメータ数は、多項式の項数にtauの分を足している。

計算結果の末尾の部分を取り出すと、 


K = 2, AIC = 64.7073005000188, AICc = 68.1358719285902, BIC = 65.9009863184139
K = 3, AIC = 66.9230544169621, AICc = 73.5897210836288, BIC = 68.5146355081556
K = 4, AIC = 52.6078217418123, AICc = 64.6078217418123, BIC = 54.5972981058042
K = 5, AIC = 54.3946154061608, AICc = 75.3946154061608, BIC = 56.781987042951
K = 6, AIC = 56.2295868715862, AICc = 93.5629202049195, BIC = 59.0148537811748

といった感じ。データは3次関数をもとにしているので、K=4が正解なのだが、各指標とも、そこが最小になっている。

該当するパラメータの推定値は、 


K = 4
# var, mean, sd, naive_se, time_series_se, geweke
alpha[1], 0.48221, 1.06165, 0.0150141, 0.028524, -0.47792
alpha[2], 0.27981, 0.55112, 0.0077941, 0.027061, 0.52308
alpha[3], -0.04185, 0.07932, 0.0011217, 0.002133, 0.99723
alpha[4], 0.12867, 0.02860, 0.0004045, 0.001369, -0.61282
tau, 0.24885, 0.13399, 0.0018950, 0.003104, -0.03633

実際は、alpha = [0, 1, 0, 0.1]、tau=0.25で作っているので、……、微妙。alpha[2]以外は、誤差の範囲内で合ってるといっていいのかな?

何はともあれ、計算できた。

[「間に合わせのパッケージ」に戻る] [「PerlとJAGSで赤池情報量基準(AIC)」に戻る] [「あえてのPerl」に戻る]