image.plot の目盛を axis 関数で描く

Rのfieldsパッケージのimage.plot関数で目盛に数値ではなくテキストを書き込みたい。axis関数を使っても何故か反映されないというバグ?がある(fieldsのバージョンは6.9.1)。例えば以下の例:

library(fields)

x <- matrix(rnorm(25),nc=5)

colnames <- paste0("col",1:5)
rownames <- paste0("raw",1:5)

# Fail
image.plot(1:5, 1:5, x, xaxt="n", yaxt="n", xlab="X-axis", ylab="Y-axis")
axis(side=2,at=1:5,labels=rownames)
axis(side=1,at=1:5,labels=colnames)

image_plot_fail

確かに最後の二行の axis 関数の結果が反映されていない(x軸, y軸の目盛に文字が書けていない)。

r- how to edit elements on x axis in image.plot
この問題はここでも議論されていて、普通のimageを使って、凡例だけimage.plotを使え、とある。

もう少し簡単に解決する方法を見つけた。

# Success
image.plot(1:5, 1:5, x, xaxt="n", yaxt="n", xlab="X-axis", ylab="Y-axis")
text(-100,-100,".")  # This is dummy, but required for drawing axis
axis(side=2,at=1:5,labels=rownames)
axis(side=1,at=1:5,labels=colnames)

image.plot と axis の間にダミーでtext関数を実行する(三行目)。以下のとおり不思議なことにこれで解決した(理由は未調査)。

image_plot_success

重点サンプリング (2)

前回は正規分布の裾の積分をなんとなく決めた提案分布による重点サンプリングで求めた。今回は提案分布の違いがどのような誤差の違いを生むのかについて実験した。ただし、今回の積分範囲は [4,∞) とした。ようするに \(f\) を標準正規分布のpdfとして

\begin{align}
\int_4^\infty f(x)dx
\end{align}

を求める問題。
“重点サンプリング (2)”の続きを読む

モデル選択の実験 (BIC を追加)

前回の記事 では AIC と AICc を比較した。今回はそれに BIC を追加してみた。BICはあまり使ったことがなかったが、個人的には結構おどろきの結果が得られた。

BIC は以下で定義される。n はデータ数、k はモデルのパラメータ数。

\begin{align}
BIC=-2\times\{\text{Maximum Likelihood}\}+(\log n)\times k
\end{align}

実際のデータ分析では当然、n は固定なので AIC とのちがいは k の前の係数が 2 という定数か、 \(\log n\) という定数か、の違いがあるが、これって同じようなもんでしょ、と思って BIC はわりとノーマークだったけど、今回実験してみて考えを改めることなった。
“モデル選択の実験 (BIC を追加)”の続きを読む

モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)

前回の「モデル選択の実験 (AIC vs LOOCV)」の続きです。

小標本の場合は、AIC じゃなくて AICc を使うといいよとのことなので、今回は、前回同様の方法で AIC と AICc を比較してみた。

真のモデルやその他のモデルの設定などは前回と全く同様。図の見方も前回と同様です。LOOCV の結果も並べてみたかったが計算量の関係で断念。

いろんな意味で手抜き気味です。あしからず。

AICc

「正規ノイズの線形モデル」のケースでは以下で定義される情報量基準。

\begin{align}
AICc = AIC + \frac{2k(k+1)}{n-k-1}
\end{align}
k : パラメータ数、n : データ数。

GLM のケースではこの定式化は使えないらしい[1]。 そんなこんなで GLM の場合は R の AICcmodavg パッケージをつかう。

glm.fit <- glm(...)
AICc(glm.fit)

とすることで AICc の値を良きに計らって計算してくれる。
“モデル選択の実験 (AIC vs AICc / R の AICcmodavg パッケージ)”の続きを読む

  1. [1] でも実は間違えてこの定式化でやってみたけどそれなりに良い結果 (AIC よりもよい結果) が得られた。

モデル選択の実験 (AIC vs LOOCV)

最近読んだ「データ解析のための統計モデリング入門――一般化線形モデル・階層ベイズモデル・MCMC (→ amazon)」に感化されて、モデル選択基準である、AIC、LOOCV (Leave One-Out Cross Validation) について実験してみた(参考:モデル選択の周辺の話の整理)。

LOOCV と AIC はある意味で両極端だ。LOOCV は利用するのに必要な仮定が少なく、直感的にも何をやっているのかわかりやすい(ただし計算量は非常に大きい)。一方で AIC は評価式の簡単さ(計算量の少なさ)とは裏腹に、背後には正則条件だとかネストだとかの直感的とは言いがたい仮定をしたりする。

まあ AIC とか LOOCV とかが常に正しいモデルを選択するわけじゃない、というのは直感的にはわかるけど、実際に試してみたことはなかったので、真のモデルを知っているという神の視点からみて、こいつらがどんな挙動を示すのかをシミュレーションで確かめてみた。

“モデル選択の実験 (AIC vs LOOCV)”の続きを読む

[R] データフレームを可視化するtabplotパッケージについて

データの列数が多くなってくる (高次元になってくる) とデータの全体像が捉えにくくなる.R-bloggers を読んでたら,Rのデータフレームを可視化するためのパッケージ tabplot なんてのがあるらしい,ということで試してみた.コレを使うとデータを効率的に可視化できるので,よくわからないデータをとりあえず可視化してみて,それからあれこれ考えてみると捗るかもしれない.

CRAN から tabplot パッケージをインストールすればすぐに使える (類似の名前で tableplot なんてのがあるので注意).
“[R] データフレームを可視化するtabplotパッケージについて”の続きを読む

Rで階段プロットを書く

知らなかったのでメモ。Rで階段プロットを書くには type=”s” を指定する。

以下は経験累積分布関数を書くためのコード片。

cumplot <- function(x,...){
  plot(sort(x),(1:length(x))/length(x),type="s",...)
}

cumplot(rnorm(100),xlab="",ylab="",xlim=c(-3,3))
curve(pnorm,add=T,col=2)
legend("topleft",legend=c("empirical CDF","normal CDF"),
       lt=1,col=1:2)

こんな図が書ける。

ほかにも type=”S” (大文字) なんてのがあって、だいたい同じなんだけど、S は上にずれる。経験累積分布の場合は小文字 s が正解。

CRPのテーブルの数の分布

Chinese Restaurant Process (→ 以前の記事) でデータ数(レストランに来る人数/壺からボールを取り出す回数)や \(\alpha\) が変化した時に利用されるテーブルの数の分布がどうなるか実験してみた(下の図をクリックで拡大)。


“CRPのテーブルの数の分布”の続きを読む

正規分布/Normal-inverse-Wishart が事後分布に収束していく様子

正規分布のパラメータ \(\mu, \Sigma\) の共役事前分布は Normal-inverse-wishart (NIW) 分布。データ数が増加した時に真のパラメータに収束していく様子を図示してみた (クリックで拡大する、かも)。


絵のせつめい。

  • データ数 n を変化させた各グラフの中で、事後分布から20組のパラメータをサンプリング。
  • 一組のパラメータを90%信頼楕円として表している
  • 点線は真の分布。黒丸は事前分布のモード(最頻値)。
  • データが少ないうちは事後分布がばらつく、つまりパラメータの不確実性が大きいが、データ数が256を超えたあたりからほぼ真の分布に収束する

下のコードのコメントにも書いたが、ハイパーパラメータを学習しない場合のレシピとして、NIW分布の平均パラメータ \(\mu_0\)、分散パラメータ \(V_0\) に自信がないときはそれぞれ \(k_0, \nu_0\) を小さく設定すればいい、はず。事前分布の影響を小さくできる。

参考その1:前回のエントリ → 逆Wishart分布を図示してみる
参考その2:事後分布のパラメータの求め方 → Wikipedia/Conjugate prior

お絵かきスクリプト in R。ellipse, MCMCpack, mvtnorm パッケージは入っていなければインストールする必要あり。

library(ellipse)   # conffidence ellipse
library(MCMCpack)  # wishart
library(mvtnorm)

#----------------------------------------------------------------------------
# Baysian estimation of 2 dimensional normal distribution
#----------------------------------------------------------------------------

d <- 2     # dimension

#----------------------------------------------------------------------------
# Hyperparameters for NIW (Normal inverse wishart)
k0 <- 0.1        # see below
mu0 <- rep(0,d)  # hyper mean
v0 <- 3.5        # see below
V0 <- diag(rep(10,d)) # hyper variance
Phi0 <- V0*(v0-d-1)
prior.hyper.par <- list(k=k0,mu=mu0,v=v0,Phi=Phi0)

#------------------------------------------------------------------------------------
# [Recipe for determining hyperparameters v0 and k0]
#  1. Less confidence for the hyper mean(mu0), take smaller k0; typically, 0<k0<<1
#  2. Less confidence for the hyper variance(V0), take smaller v0; typically, v0~1+p
#------------------------------------------------------------------------------------


#-----------------------------------------
# Bayesian estimation ( bayesian update)
bayes.update <- function(X,hyper.par){

  k0 <- hyper.par$k
  mu0 <- hyper.par$mu
  v0 <- hyper.par$v
  Phi0 <- hyper.par$Phi

  if(!is.matrix(X)){
    # This is n=1 case that we have to deal with as special because of R issue
    n <- 1
    EX <- X
    X <- t(X)
    C <- matrix(rep(0,d*d),nc=d)
  } else {
    n <- dim(X)[1]
    EX <- colMeans(X)
    C <- (n-1)*cov(X)
  }
  mu <- (k0*mu0 + n*EX)/(k0+n)
  k <- k0+n
  v <- v0+n
  Phi <- Phi0 + C + k0*n/(k0+n)*((EX-mu0) %*% t(EX-mu0))

  list(k=k,mu=mu,v=v,Phi=Phi)  # This new parameters list is a bayesian update result!!

}

#------------------------------------------------------------------
# Sampling from posterior
sample.posterior <- function(n,bayes.fit){
  k <- bayes.fit$k
  v <- bayes.fit$v
  Phi <- bayes.fit$Phi
  mu <- bayes.fit$mu
  result <- list()
  for(i in 1:n){
    # Sampling mu and V from NIW
    V <- riwish(v,Phi)  # 1. sampling covariance matrix by inverse wishart
    result[[i]] <- list()
    result[[i]]$V <- V
    result[[i]]$mu <- rmvnorm(1,mu,V/k) # 2. sampling mu by normal dist.
  }
  result
}

#--------------------------------------------------------------
# Test case
#  --- true distribution
V.true <- matrix(c(3,-1.5,-1.5,1.8),nc=d)
mu.true <- rep(10,d)
#  --- generating random variables
X <- rmvnorm(20000,mu.true,V.true)


par(mfrow=c(3,4))  # dividing graphic device
j <- 1
ids <- 2^seq(0,11)
hc <- rainbow(12)
for(i in ids){
  Y <- X[1:i,]
  bayes.fit <- bayes.update(Y,prior.hyper.par)  # estimating posterior
  post.samp <- sample.posterior(20,bayes.fit)  # sampling from posterior
  plot(Y,pch=".",col="gray",xlim=c(-10,20),ylim=c(-10,20),xlab="",ylab="") # just plotting raw data

  # plotting posterior samples
  for( theta in post.samp ){
    mu <- theta$mu
    V <- theta$V
    elp <- t(apply(ellipse(V,level=0.9),1,function(x) x+mu))
    lines(elp,col=hc[j],lw=2)
  }

  # plotting true distribution
  true.dist <- t(apply(ellipse(V.true,level=0.9),1,function(x) x+mu.true))
  lines(true.dist,type="l",lw=2,lt=4)

  # plotting mode of the prior
  prior <- ellipse(Phi0/(v0+3), level=0.9) + mu0
  lines(prior,lw=2)

  text(20,-7,paste("n =",i),pos=2,cex=2)
  j <- j+1
}

dev2bitmap(file="normal_posterior.jpg",taa=4,gaa=4,width=15,height=12)

逆Wishart分布を図示してみる

正規分布のベイズ推定で共分散行列の(共役)事前分布として使われる逆Wishart分布(→ wikipedia:en)というものがありますが、あまりイメージしにくかったので図にしてパラメータの変化がどのような違いをもたらすか調べてみました。

といっても逆Wishart分布は n x n 行列を値にとる確率分布なので、2×2 行列を考えても3次元なので(※対称行列)そのままプロットするとわかりにくい。そこで ellipse パッケージを使って可視化してみます。
“逆Wishart分布を図示してみる”の続きを読む