モデル選択の周辺の話の整理

モデルを選択したり、変数を選択したり、というようなことに関係しそうなネタを簡単に整理してみた。

情報量基準

AIC / BIC / DIC / TIC のような。データのあてはまりのよさとモデルの複雑度を天秤に図るタイプのやつ。

たくさん種類があるのは確率モデルに関する仮定と汎化誤差の近似の仕方の違いによるものだと理解している。

検定

回帰係数が0である、という帰無仮説を検定することである変数が貢献しているかどうかを定量化するタイプのやつ。棄却されなければ「えい!」と変数を削ってしまう。

ささっと分析してデータの雰囲気を掴みたい時に使うことはある。

L1正則化

寄与度の小さな(ある閾値より小さな)係数をゼロにしてしまう、というような感じのやつ。

事前分布としてラプラス分布を使うことに相当。単純に寄与度が低い変数は消してしまえ!というノリなのだろうか?もっと深遠な背景があるのだろうか?勉強不足でよくわからない。

ベイズモデル選択

複数のモデルに事前分布を設定して、「モデルの事後分布」を計算するたぐいのもの。事後分布が求まったあとはMAPなものを選んでくるか、事後分布で平均をとってしまうか。

たとえばディリクレ混合過程。これはGMMのような混合モデルの混合数の事後分布を求めることができる。

PRMLの変分ベイズのところで出てきた関連度自動決定もこのタイプだと思っていいのだろうか。これも勉強不足により不明。

***

いろいろと抜けがあるとは思うが、とりあえずすぐに思いついたのはこれくらい。場合によっては追記します。

ハイブリッドモンテカルロの実験

相関の強い二変量正規分布に対してハイブリッドモンテカルロを使ってみた。上から順に、サンプリング結果、x1の自己相関、x2の自己相関。

自己相関ほぼ完全になし、という結果になった。ギブスサンプラーだとこうはいかない。ただし、

  • 微分方程式を解く時間のスケールが小さすぎると自己相関が出たので良い感じのスケールをちょっとだけ探索した。
  • ステップ幅を固定にしたら怪しげな自己相関の挙動がでた。
  • 計算時間でギブスサンプラーと比較してどちらが有利かは今回は検討してません。

という点は追記しておきます。

“ハイブリッドモンテカルロの実験”の続きを読む

対数線形モデルとエントロピー最大化の関係

昔の勉強ノートを引っ張りだしてくるシリーズ.

機械学習の対数線形モデルが最大エントロピー法とも呼ばれる,みたいな記述は頻繁に目にするし,統計力学のボルツマン分布の話とか考慮すれば,なんとなくそうなってそうな気はするけど,実際どうなの?というのを (たんに好奇心を満たすために) 調べてみた.実用上は何の意味ないと思う.

対数尤度関数に L1 正則化項を加えるタイプの目的関数を使った場合,もはやエントロピーは最大化されない,とかそういうわりとどうでもいいことがわかったりするかもしれない.

概要

「言語処理のための機械学習入門 (→ amazon) 」などに出てくるタイプの対数線形モデルの係数の最尤推定量が,エントロピーを「ある制約条件下」で最大化した場合のラグランジュ未定乗数に対応することを説明する (クロス表の対数線形モデルとはたぶん別物).

ただし,記号が煩雑になるのを避けるため,対数線形モデルとほぼ同一の構造を持ち,記号が煩雑でない条件付きロジットモデルがエントロピー最大化と等価であることを見る.

本文の最後に対数線形モデルと等価なエントロピー最大化問題を示す.多少ややこしくなるが,同じ方針で証明可能.
“対数線形モデルとエントロピー最大化の関係”の続きを読む

チューリングマシンと限定合理性 : 「行動ゲーム理論入門」を読んだ

この本「行動ゲーム理論入門」はたまたま本屋で見かけてパラパラと見ていたら、経済学の本にもかかわらず「チューリングマシン」だとか「強化学習」だとかいう一見経済学とは関連の薄そうな単語があったので、興味深いな、と思って脊髄反射的に購入した。僕はこの分野は全く知らない状態でこの本を読み始めたのだけど、非常に刺激的な本だったので記憶が鮮明なうちに書いておくことにする。

“チューリングマシンと限定合理性 : 「行動ゲーム理論入門」を読んだ”の続きを読む

Google検索時に出てくるevernoteのアレの位置を変更する

Chrome で evernote プラグインを使ってるときに google 検索すると evernote の検索結果も同時に出てくるアレの話です。ピンとこない人は関係ないはずです。

あの検索結果は個人的には結構便利だとは思うのですが、レイテンシがわりと大きいため、一番上のリンクをクリックしようとした瞬間にあの水色のボックスが出現したりすることが多いのが少し困る。

かといって機能を切ってしまうのもなんとなく忍びないのですこし改良してみた。最初は javascript だけでやろうとしたけど、手持ちの chrome には stylebot という拡張機能がついていたのでそれを使うことにした。
“Google検索時に出てくるevernoteのアレの位置を変更する”の続きを読む

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のテーブルの数の分布”の続きを読む

Rにおける値渡しと参照渡し (2)

昨日の記事を書いたあと、twitter で tracemem 使うといいよ、と教えていただきました (@sfchaosさん、ありがとうございました!)。

この関数ははじめて知ったのですが、ヘルプを見て意訳すると以下の様な感じです。

tracemem(x) を実行すると x について duplicate (複製)が生じた時にメッセージを表示する。これは2つのオブジェクトがメモリを共有している場合において、1つが変更された時に生じる。ちなみに untracemem(x) でメッセージフックを解除できる。

ということで関数の内部でコピーが起きたか、を判別するのにピッタリです。

今回の内容:

  • tracemem を使って値渡し、参照渡しを確かめる
  • 参照渡しなのか、値渡しなのか、が確定するタイミングについて考える

tracemem を使った方法

以下の方法は @sfchaos さんに教えてもらった方法そのままです。前回同様、prod1 は行列 A について値渡しになりそうな関数で、prod2 は参照渡しになりそうな関数です。

prod1 <- function(A,x){
  A <- A + diag(x)
  A %*% x
}

prod2 <- function(A,x){
  A %*% x
}

N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

tracemem(A)
# => [1] "<0x7e790008>"
invisible(prod1(A,x))  # invisible は結果を print しないようにする関数
# => tracemem[0x7e790008 -> 0x7dfe0008]: prod1  # コピーが発生した!

invisible(prod2(A,x))
# => 何も表示されない!=コピーが生じていない!

ということで、前回は実行時間から推論しただけでしたが、やはり引数を変更するとコピーが生じる、ということで間違いないことが確認できました。

複製はどのタイミングで生じるのか?

前回、以下のように書きました。

そして、引数が変更されるかされないかはパースした段階でわかる(なのでパースの段階で値渡しか参照渡しかを判別することが可能)

ですが、これは誤りでした。実際、以下のような関数はパースの段階で値渡しか参照渡しかを判別できるでしょうか?

prod3 <- function(A,x,add=T){
  if(add) A <- A + diag(x)
  A %*% x
}

add が真のときは引数が変更されて、偽のときは引数が変更されません。

これをパースの時点で判別しようとすると、Rの副作用を作らないという原則から add がいかなる値であろうともコピーを行うことになるはずです。

もう一つの可能性としては、実際に変更が起きたその瞬間にコピーを作る、というものがありえるでしょう。

実験してみます。

N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

tracemem(A)
# => [1] "<0x7dfe0008>"

invisible(prod3(A,x,add=T))
# => tracemem[0x7dfe0008 -> 0x7ef40008]: prod3

invisible(prod3(A,x,add=F))
# => 何も表示されない!=コピーが生じていない!

同じ関数でも引数の状態によってコピーが発生したり、しなかったり。add=F のときは prod3 の if の内部まで進まないため、A が変更されず、したがってコピーが発動しない、ということになります。

つまりRは

  • 引数のコピーを作るか(値渡しか)、作らないか(参照渡しか)は実際に引数が変更される瞬間ギリギリまで判別しない、
  • 変更される瞬間(直前?)で重い腰を上げてコピーを作成する(遅延評価)

という挙動をしているようですね。

まとめ

  • Rは基本的には値渡しであり、C++のように引数として与えられた変数を変更することで外側の世界に影響させることはできない
  • 関数内部で引数を変更しない場合は、コピーが生じない (C++ の const& のようなイメージ)
  • 関数内部に引数を変更するコードがあったとしても、実際に引数が変更される段まではコピーは生じない (遅延評価)

Rにおける値渡しと参照渡し

Rの関数に引数を渡すと値渡しになる、とずっと信じていたわけだけど、どうも違うらしい。Rはどうやら「自動的に」値渡しすべきか、参照渡しにすべきか、を判断しているようだ。

C++みたいな言語では「フツーに」引数を渡すと全部値渡しになって(特に行列のような巨大なオブジェクトを渡す場合には)効率が悪いので、ポインタや参照で渡したり、副作用を気にする場合は const 参照で渡したりする。

さて、上で述べた R の自動判断機能は以下のような事実に基づくようだ。

  • そもそも値渡しは関数に引数を渡した段階でオブジェクトのコピーを生成して「引数として渡された変数を関数の内部で変更してもスコープの外では値が変更されない」ことを保証するのだが、
  • オブジェクトが関数の内部で変更されないことが保証されるならば「関数の実行中はスコープの外でも値が変更されない」ことは保証されるのでコピーをそもそも生成する必要がない、
  • そして、引数が変更されるかされないかはパースした段階でわかる(なのでパースの段階で値渡しか参照渡しかを判別することが可能) → 間違いでした。詳細はRにおける値渡しと参照渡し(2)に書きました。

実験

これは以下のようにして確かめることができる、はず。

prod1 <- function(A,x){
  A[1,1] <- A[1,1] + 1
  A %*% x
}

prod2 <- function(A,x){
  x[1] <- x[1] + 1
  A %*% x
}
  • prod1, prod2 はともに引数として渡された行列 A とベクトル x の積を計算する
  • prod1 では A[1,1] に 1 を加える
  • prod2 では x[1] に 1 を加える(これはprod1の代入作業のコストと揃えるため)
  • その後、積を計算する

ということをやっているのだが、上で述べたようなことがただしければ、

  • prod1 では行列 A のコピーが発生し
  • prod2 ではベクトル x のコピーが発生する

ため、より巨大なオブジェクトを渡すことになる prod1 が (生じる四則演算の数は同一にもかかわらず) 大幅に遅くなることが予測される。

実際、1000 x 1000 程度の行列を考えると以下の様な結果が得られる。


N <- 1000
A <- matrix(rnorm(N*N),nc=N)
x <- rnorm(N)

system.time( for(i in 1:1000) prod1(A,x) )
# =>   user  system elapsed
# =>  11.03    3.50   14.56

system.time( for(i in 1:1000) prod2(A,x) )
# =>  user  system elapsed
# =>  4.59    0.01    4.62

prod1のほうがかなり遅くなっていることが分かる。これより、上で述べたことはおそらく正しいだろうと結論される。

まとめ

このことを知ったからといって高速なプログラムが書けるようになるわけではないが、少なくとも「こんな大きなオブジェクトを関数に渡すのは気が引ける(なのでグローバル変数にしてしまえ)」という不安を払拭することができると思う。

【追記】続きを書きました。

正規分布/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)