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分布を図示してみる”の続きを読む

指数型分布族のメモ

モチベーション

推定や検定の理論を扱うときに、正規分布の場合~、ポアソン分布の場合~、みたいに分布ごとに性質を調べてたんじゃやってられない、ということで、さまざまな密度関数をある種のクラスとしてまとめて考えることで、一般論を展開したい。そのための枠組みの一つが指数型分布族。

実際、指数型分布族に該当する分布は推定や検定などで「極めていいかんじの」性質を持っている。一方で実用上よく使う分布の多くが、指数型分布族というクラスを「継承」している(ほとんどが、と言ってしまってもいいかもしれないが、機械学習系の確率モデルには当てはまらないことが多い)。

要するに、指数型分布族の理論について知っていれば、実用上よく使う分布たちが極めてよい性質を持っていることが統一的に理解できる。ちなみに、GLM(一般化線形モデル)も基本的には指数型分布族前提で定式化されている。

適当に書いて見ましたが、こんなところでしょうか。。。

定義

n次元の確率変数 \(X\) が d-パラメータの指数型分布族に従うとは、あるパラメトリゼーション \(\boldsymbol{\theta}\) と d 個の関数 \(T_i:\mathbb{R^n}\to\mathbb{R}\hspace{5pt}(i=1\cdots d)\) が存在して、\(X\) の密度関数が
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)
\end{align*}
と書けること。\(h, A\) はそれぞれスカラー値をとる関数。\(T_i(\cdot)\) は十分統計量と呼ばれる。

指数型分布族の分布の例

正規分布、ポアソン分布、指数分布、ガンマ分布、ディリクレ分布、負の二項分布などなど。数え出したらきりがない。

指数型分布族でない分布の例

一様分布、コーシー分布など。混合ガウス分布とか隠れマルコフモデルとか、その手のいわゆる「特異モデル」は該当しない。

Natural exponential family

分布族を
\begin{align*}
f(x|\boldsymbol{\theta}):=h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot x-A(\boldsymbol\theta)\big)
\end{align*}
に限定したものが natural exponential family。

平均と共分散の計算

指数分布族に従う d 次元確率変数 \(T(X)=(T_1(X),\cdots,T_d(X))\) の平均と共分散は
\begin{align*}
E(T_i(X))&=\frac{\partial A}{\partial\theta_i}\\
Cov(T_i(X),T_j(X))&=\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
となる。なぜなら、正規化条件
\begin{align*}
1=\int h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
\end{align*}
の両辺を \(\theta_i\) で微分すれば
\begin{align*}
0=E_\theta(T_i(X))-\frac{\partial A}{\partial\theta_i}\cdot1
\end{align*}
となり、もう一度、今度は \(\theta_j\) で微分すれば
\begin{align*}
0=&\int T_i(x)T_j(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx\\
&-\frac{\partial A}{\partial \theta_j}\int T_i(x)\cdot h(x)\cdot\exp\big(\boldsymbol{\theta}\cdot\boldsymbol{T}(x)-A(\boldsymbol\theta)\big)dx
-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&E(T_i(X)T_j(X))-E(T_i(X))E(T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}\\
=&Cov(T_i(X),T_j(X))-\frac{\partial^2 A}{\partial\theta_i\partial\theta_j}
\end{align*}
ということでOK(積分と極限の交換の条件とかは無視してしまいましたが)。

Chinese Restaurant Process を発生させてみた

タイトルの通りです。単なる遊びです。Blackwell-McQueen Polya’s urn scheme とか表記されたりもするみたいです。

Chinese Restaurant Process (CRP)

よく(名前にもあるように)中国料理屋の例が出ますが、個人的には「初期の壺に連続的なスペクトルの色(虹色)を持つボールが入っているバージョンのポリアの壺」と考えるほうがわかりやすかったです。つまり、以前に記事にしたポリアの壺の色数が(可算とは限らない)無限集合に拡張されたもの、ということです。兎に角、CRP では以下のスキームで無限の長さの確率変数列を発生させます。

  1. \(G_0\) を適当な確率分布、\(\alpha\) を正の実数とします
  2. \(k\leftarrow1\)とする
  3. \(G_0\) から乱数(色)をひとつサンプリング
  4. その乱数を記録し、同じ色のボールを空の壺にいれる
  5. \(k\leftarrow k+1\)とする
  6. \(\alpha/(\alpha+k-1)\) の確率で \(G_0\) からサンプリングし、それ以外のときは壺からボールを一つ取り出してそれをサンプルとする
  7. サンプルと同じ色のボールを壺に入れる
  8. ステップ5に戻る

実際には無限に数列をつくることはできないので、10000個サンプルした段階でやめることにします。長さ10000の系列を15回生成して、その分布をプロットしてみます。\(G_0\) は標準正規分布、\(\alpha=3\) としました。

# Generating n random variables via CRP(alpha, G0)
rcrp<-function(n,alpha,G0){
  res<-numeric(n); res[1]<-G0(1);
  for(i in 2:n){
    p<-alpha/(alpha+i-1)
    res[i]<-ifelse(runif(1)<p,G0(1),sample(res[1:(i-1)],1))
  }
  class(res)<-"crp"
  res
}

# Simple simulation
par(mfrow=c(5,3)); par(mar=rep(1.5,4)); par(oma=c(1,1,1,1))

N <- 10000

for(i in 1:15){
  plot(table(rcrp(N,3,rnorm))/N,xlim=c(-3,3),ylim=c(0,0.4),ylab="",xaxt="n")
  axis(side=1,at=-2:2)
  curve(dnorm,add=T)
}

ここからわかること:

  • \(G_0\) は連続値をとる密度関数だが、得られる分布は離散的(実はこれは10000個で生成をやめずに、無限個生成しても離散的)
  • 10000個のサンプルのほとんどは数~十数のポイントに落ちる

いろいろ調べると(証明すると)わかること:

  • 実は得られた15個の離散分布はディリクレ過程(←ここでは説明していない)に従う「ランダム確率分布」である
  • (これは、もちろん有限で止めているという意味では正確な表現ではない)
  • 本当のディリクレ過程からのサンプルされる分布は(離散だけども)サポートは無限集合になる
  • CRPによって得られる確率変数列は過去に依存しているので、iidには見えない。しかし、実際にはこれらは得られた離散分布からのiid とみなすことができる

ディリクレ分布に関するメモ (6)

ディリクレ分布に関するメモ (3) ではガンマ乱数を正規化することでディリクレ分布に従う乱数を発生させました。

今回は別の方法を試してみます。

ポリアの壺

ポリアの壺(Polya’s urn)というのは以下のような手続きのことです。

  1. n種類の色のボールが入った壺を考える
  2. 壺から一つボールを取り出す
  3. 取り出したボールを壺に戻し、同じ色のボールを壺に追加する
  4. ステップ2に戻る

この手続きによって、壺の中のボールは増加し続けます。その色のシェアはランダムに変動しますが、やがて収束します。この収束した時のシェアを \(p^*\) とします。この極限のシェア \(p^*\) はボールを取り出した経路に依存し、これ自体が確率変数です。これは、例えば壺に赤と青のボールが一つずつ入っていたとして、最初の数回でたまたま赤が連続して出ると、その後も赤が出やすくなり、極限では赤のシェアが90%ということも起こりえるわけです。

さて、n色のポリアの壺を考えます。最初に \(\alpha=(\alpha_1,\cdots,\alpha_n)\) 個のボールが壺に入っていたとすると、このポリアの壺の極限分布が \(\text{Dirichlet}(\alpha)\) に従います(ちなみに、これは \(\alpha\) が自然数でなくても成立します)。
“ディリクレ分布に関するメモ (6)”の続きを読む

ディリクレ分布に関するメモ (5)

以前、ディリクレ分布の平均、分散、共分散を計算しました。

ところが今日、ディリクレ過程に従うランダム確率分布がなぜ離散的か?という問題に対する証明を与えている

Ferguson Distributions Via Polya Urn Schemes; David Blackwell and James B. MacQueen, Ann. Statist. Volume 1, Number 2 (1973), 353-355.

という論文を読んでいて任意の高次モーメントが必要になったので、計算してみました。
“ディリクレ分布に関するメモ (5)”の続きを読む

ディリクレ分布に関するメモ (4)

前回はディリクレ分布にしたがう乱数を発生させてそれを可視化しましたが、乱数発生のアルゴリズムがなぜ正当化されるのかについては説明しませんでした。

今回はそのあたりについて書いてみます。前回書いたように、パラメータ \(\alpha\) をもつディリクレ分布からのサンプリング方法は以下のようになります。

  1. \(i=1,\cdots,d\) に対してシェープパラメータが \(\alpha_i\) 、スケールパラメータが1のガンマ分布を発生させる。これを \(X_i\) とする。
  2. \(X=X_1+\cdots+X_d\) として、\(Y_i=X_i/X\) としてこの \(Y\) を記録する

このように定義した \(Y_i\) がディリクレ分布に従うことを以下で証明します。 “ディリクレ分布に関するメモ (4)”の続きを読む

ディリクレ分布に関するメモ (3)

すこし予定を変更して、今回はディリクレ分布を可視化してみます。
次元が高いとプロットできないので、3次元のディリクレ分布を可視化することにします。


三次元なのでこのような三角形の上に確率変数が分布します。ディリクレ分布のパラメータをいろいろ変えながら、この三角形上にランダムにサンプルされたディリクレ分布をプロットしていきます。そのためにはディリクレ分布に従う乱数を発生させるアルゴリズムが必要になります。 “ディリクレ分布に関するメモ (3)”の続きを読む

ディリクレ分布に関するメモ (2)

前の記事 ではディリクレ分布の正規化定数を計算しました。

ディリクレ分布は正規化定数を含めて以下のようになりました。
\begin{align*}
&f(x)=\frac1{B(\alpha)}\prod_{i=1}^d x_i^{\alpha_i-1}\qquad x\in \Delta\\
&\text{where} \qquad B(\alpha)=\frac{\prod_{i=1}^d\Gamma(\alpha_i)}{\Gamma(A)}
\end{align*}

今回は平均、分散、共分散を計算します。これらの計算には正規化定数の計算と類似点が非常に多いため、前回の記事を読んでない方は先にそちらをどうぞ。 “ディリクレ分布に関するメモ (2)”の続きを読む