Rで将来のPVをMCMC推定する

学習中のベイズ推定。
特にMCMC推定はなんだかしっくりくるので本を読んでみると、難解な数式やらなにやらと格闘。
なにはともあれちょっとやってみたい!ということでひとまずやってみました。

今回は色々あるベイズ推定関連パッケージからMCMCpackを選択。

ブログ「Taglibro de H」の「MCMCの勉強(1) [統計]」を参考にさせていただきながらやってみました。

データは私のブログのアクセスデータを利用。
アクセスデータの中からPVを利用し、今後のPVを推定していきます。

まずは、パッケージの読み出し。
> library(MCMCpack)

初期化。
> set.seed(1)

平均10、標準偏差3の乱数を1000個生成してXに代入。
> m <- 10 > s <- 3 > x <- rnorm(1000, m, s) 事前分布を作成していきます。 MCMC推定に使用する関数を用意。 > llnormfun <- function(beta, x) { if (beta[2] < 0.0) l <- -Inf else l <- sum(log(dnorm(x, mean=beta[1], sd=beta[2]))) return(l) } 事後分布を作成していきますが、ここでブログのPVデータから平均と標準偏差を算出します。 PVを変数dataAに格納しています。 mean()関数で平均を、sd()関数で標準偏差を出していきます。 ※下記ソースではA.PVという変数を別途作成していますが、これは元データにPV以外が入っているためです。特に意味はないです。 > A.PV <- dataA > A.mean <- mean(A.PV) > A.mean
[1] 41.14953
> A.sd <- sd(A.PV) > A.sd
[1] 14.99736

上記データより、推定値の初期値として平均を変数「A.mean」、標準偏差を変数「A.sd」とし、繰り返し10000回、初期値依存として捨てる数を500として、MCMC推定を実行します。

> post.samp <- MCMCmetrop1R(llnormfun, theta.init=c(A.mean,A.sd), x=x, thin=1, mcmc=10000, burnin=500, tune=0.1, verbose=500, logfun=TRUE) 上記を実行すると、長々と詳細データが表示されます。 さっそくプロット。 > plot(post.samp)

スクリーンショット 2013-06-13 1.55.37

summary()関数でデータの概要をチェック。

> summary(post.samp)

Iterations = 501:10500
Thinning interval = 1
Number of chains = 1
Sample size per chain = 10000

1. Empirical mean and standard deviation for each variable,
plus standard error of the mean:

Mean SD Naive SE Time-series SE
[1,] 31.17 5.079 0.05079 4.1754
[2,] 18.91 0.963 0.00963 0.5957

2. Quantiles for each variable:

2.5% 25% 50% 75% 97.5%
var1 23.03 26.58 30.79 35.81 39.73
var2 16.20 18.74 19.18 19.56 19.97

この推定が正しいか否かを確認してみます。
実は確認用にアクセスデータを作る際に、直近の220日分を抽出し、前半と後半を分けていました。
今回仕様したはのは前半つまり昔のデータを使ったので、そのまま将来の数字である後半部と比較してみます。

後半データの平均と標準偏差を出してみます。
後半データを変数dataB(※2列目がPVデータです)としています。

> mean(dataB[,2])
[1] 29.55455

> sd(dataB[,2])
[1] 15.54844

おお!
なかなか良い感じになりました。

しかしながら、数値の見方がわからないところが多々あるので本当にハマっているか否かは自信ないなぁ。理論をもう少し学んで実戦レベルにもっていきます。

Leave a Reply

%d人のブロガーが「いいね」をつけました。