2019-09-24

特異モデルに対してBridge Samplingで周辺尤度の推定は可能か

暫定的な見解：難しいのでは？（[編集後記 2019/9/28] warp-U bridge samplingを使うとうまくいくかもしれない。推定の分散が他の方法と比べてどれだけ小さいかが気になるところ）

定義：周辺尤度

まず言葉の定義をすると、モデルを $p(X|\theta)$ , 事前分布を $p(\theta)$ としたとき、周辺尤度 $p(X)$ を
$p(X):= \int p(X|\theta) p(\theta) d\theta$
と定義する。

周辺尤度の評価方法

周辺尤度は積分評価なので以下の評価が考えられる：

積分を陽に評価
数値積分(sampling technique)で評価
積分を漸近的に評価

1は古典的なアプローチでは難しいのはみなさんご存じの通り。現代的なアプローチとして代数幾何を使うものもある[1]が一般的に適用するには難しい。
2は次元が高い時には次元の呪いで誤差が累積したり分散が大きくなってしまったりします。
3はBIC, WBIC, sBICなど情報量規準的なアプローチ。

この記事の目的

この記事は上記の2の評価の1つであるbridge samplingが高い評価を得ているようなので簡単な検証をしてみました。
bridge samplingについては岡田先生のこちらのreviewを参考にするとよいかと思います。具体的な計算方法は清水先生のこちらの記事を見てください。どうやら心理学系の方でよく使われているようです。この記事を書こうと思ったモチベーションの1つは清水先生のこちらのスライドの中でWBICとbridge samplingではおそらくbridge samplingの方が精度が良いであろうというところに関する検証になります。

正則モデルのbridge sampling

まずexactな周辺尤度が知られている以下のケースでbridge samplingが正しく推定できるかを見てみます。radiata pine datasetというデータを使って以下の回帰モデルを考えます:
$y_i= \alpha + \beta (x_i - \bar{x}) + \epsilon_i, \quad \epsilon_i \sim N(0, \tau^{-1})$
事前分布の設定についてはこちらの論文に従いました。exactな対数周辺尤度は論文内にも記載があり、また著者のHPにもその計算プログラムが公開されています。
bridge samplingは以下で計算させました。

#radiata.stan
data {
  int<lower=1> N;
  vector[N] Y;
  vector[N] W;
  real meanW;
}

parameters {
  real<lower=0> tau;
  real alpha;
  real beta;
}

model {
  target += gamma_lpdf(tau | 6, 4*300*300);
  target += normal_lpdf(alpha | 3000, 1.0/(0.06*sqrt(tau)));
  target += normal_lpdf(beta | 185, 1.0/(6*sqrt(tau)));
  for(i in 1:N)
    target += normal_lpdf(Y[i] | alpha+beta*(W[i]-meanW), 1.0/sqrt(tau));
}

library(rstan)
library(bridgesampling)

lm.model <- stan_model(file='radiata.stan')
d <- read.table("RadiataPine.csv", sep=",", header=T)
data_bs <- list(N=nrow(d), Y=d$y, W=d$x, meanW=mean(d$x))

fit <- sampling(lm.model, data=data_bs, iter=10000, chains=4)
bs <- bridge_sampler(fit, silent = T)

結果は以下の通りでした：

手法	評価値
exact	-310.128
bridge sampling	-310.679

少しずれていますがこんなもんでしょうか。

特異モデルのbridge sampling

次に特異モデルでbridge samplingをしてみましょう。ここでは混合数2の混合正規分布を考えます:
$\alpha N(\mu_1, 1) + (1-\alpha) N(\mu_2, 1)$
データ生成分布は $N(0, 1)$ とします。
このケースだとexactな周辺尤度は知られていない（はずな）ので、Nealの方法＆プログラム[2]で周辺尤度を評価してbridge samplingの評価値との比較を行います。

# normalmix.stan
data {
  int<lower=1> N;
  vector[N] Y;
}

parameters {
  real<lower=0, upper=1> a;
  ordered[2] mu;
}

model {
  mu[1] ~ normal(0,sqrt(10));
  mu[2] ~ normal(0,sqrt(10));
  for(n in 1:N)
    target += log_sum_exp(
      log(1-a) + normal_lpdf(Y[n] | mu[1], 1),
      log(a) + normal_lpdf(Y[n] | mu[2], 1)
    );
}

library(rstan)
library(bridgesampling)

N <- 50
Y <- rnorm(N, 0, 1)
gmm.model <- stan_model(file='normalmix.stan')

data_1 <- list(N=N, Y=Y)
fit_1 <- sampling(gmm.model, data=data_1, iter=10000, chains=3)
BR_tmp <- bridge_sampler(fit_1, silent = T)
bridge_res <- BR_tmp[[1]]

上記のコードをベースに100回のiteratoinでNealのMonte Carlo evaluation, bridge sampling, WBIC, WBICを補正した値, sBICの平均で比較してみたのが以下の結果です:

手法	mean
Monte Carlo evaluation	-74.11
bridge sampling	-70.70
WBIC	-73.18
WBIC - $\hat{\nu}(t_w)$	-73.72
sBIC	-75.61

このケースであればNealのtechnical reportと同じでパラメータは3次元であるので最初に述べた次元の呪いの影響はあまりないと考えられるためMonte Carlo evaluationが最も精度が高いと考えてよいと思います。それと比べるとbridge samplingは結構overestimationしているように見えます。WBICもoverestimationすると指摘されていますが、それより大きくなっています。

特異モデルのbridge samplingはなぜうまくいかないか

（コードが正しいとするならば）実装で提案分布がnormal distributionになっていることが原因ではないか（要検証）。bridge samplingの提案分布は"事後分布と類似していて重なりが大きいように選ばれる分布"[3]ということだが、特異モデルの事後分布は事前にどうなるか予想が難しいと思うので（というのもデータ生成分布に依存してしまうので）、良い感じの提案分布を一般に用意するのは難しいと思う。となると、特異モデルでbridge samplingをうまく使うのは難しいのではないか。

bridge samplingなにもわからんマンなので、間違っていたらコメントをしていただけるとありがたいです。

[編集後記 2019/9/28] [1609.07690] Warp Bridge Sampling: The Next Generation なら特異モデルでもうまくいくかもしれない。まだpreprintだからか実装コードが公開されていないので公開されたら試してみます。

参考文献：
[1] LIN, S. STURMFELS, B. AND XU, Z. (2009). Marginal likelihood integrals for mixtures of independence models, JMLR (10), 1611–1631.
[2] NEAL, R. (1999). Erroneous results in ‘Marginal likelihood from the Gibbs output’. Tech. rep., Unpublished. University of Toronto, Toronto.
[3] 岡田謙介 (2018). ベイズファクターによる心理学的仮説・モデルの評価. Japanese Psychological Review, Vol. 61, No. 1, 101–115.
[4] 階層ベイズと自由エネルギー. http://norimune.net/3107

2018-06-12

実対数閾値の定義に使われているゼータ関数はなぜゼータ関数と呼ばれているのか

実対数閾値の発表をするとたまにタイトルのことを聞かれることがあるのでそれについてのメモです。

モデルの分布を $\displaystyle p(x|\theta)$ 、真の分布を $\displaystyle q(x)$ とした時、その間のKullback-Leibler divergence $\displaystyle K(\theta)$ は
$\displaystyle K(\theta) := \int q(x) \frac{q(x)}{p(x|\theta)} dx$
と定義されます。

ここで特異学習理論においてゼータ関数は、事前分布 $\displaystyle \varphi(\theta)$ 、 $\displaystyle z \in \mathbb{C}$ に対して、
$\displaystyle \zeta(z) := \int K(\theta)^z \varphi(\theta) d\theta$
と定義されます。

このゼータ関数の最も0に近い極の符号を反転させたものを実対数閾値 (real log canonical threshold)や学習係数 (learning coefficient)などと呼ばれていて、特異学習理論においては重要な役割を果たしています。

さて、上記の関数はなぜ「ゼータ関数」と呼ばれているのでしょうか。我々のよく知っているリーマンゼータ関数は
$\displaystyle \zeta(s) := \sum_{n=1}^{\infty} \frac{1}{n^s}$
という形をしていて、全く違うように見えます。

もちろんリーマンゼータ関数は色々な表現方法があり、たとえば
$\displaystyle \zeta(s) = \int_{0}^{\infty} x^s \frac{1}{\Gamma(s)x(e^x - 1)} dx$
と書けば、少し近い形にはなりますがこの表現が元でゼータ関数と名付けられたとは考えにくい。

話がいきなり飛んで答えに近づきますが、1954年にGelfandが次のような予想をしました：
For $\displaystyle s \in \mathbb{C}$ and $\displaystyle f \in \mathbb{R}[x_1,...,x_n]$ , given a test function $\displaystyle \Psi$ , the analytic map on $\displaystyle \{s \in \mathbb{C} | Re(s)>0 \}$ ,
$\displaystyle s \mapsto \int_{\mathbb{R}^n} |f(x)|^s \Psi(x) dx$
has a meromorphic continuation to $\displaystyle \mathbb{C}$ .

この問題に対して、AtiyahそしてBernstein-Gelfandがそれぞれ特異点解消定理を用いて証明を与えました。さらに後にBernsteinが前回記事でも名前を触れたD加群によって、特異点解消定理を用いない証明を与えました。

また、これと似たものでIgusaによりp進数版のゼータ関数：
$\displaystyle Z_f(s) := \int_{\mathbb{Z}_p^n} |f(x)|_p^s dx$
も作られ、ある方程式のmod $\displaystyle p^k$ の解の数の母関数と見ることができるようだ。詳しくはIgusa先生の解説記事[1]が日本語で無料で手に入るのでそちらを読みましょう。

今日このタイプの関数は局所ゼータ関数、井草ゼータ関数、Gelfandゼータ関数と呼ばれており、学習理論で使われているものもこのタイプのゼータ関数です。

参考文献
[1] 井草凖一(1994), 局所ゼータ関数について, 数学 vol.46, No.1 ,23-38.

2018-06-11

ホロノミック勾配法勉強メモ

ホロノミック勾配法の勉強メモです。

ホロノミック勾配法のよく知られている応用として

・確率分布の正規化定数の評価

・何らかの最適化

がありますが、本記事は前者についての個人メモです。後者についてはデンソーITラボの
ホロノミック勾配法 – Tech.D-ITlab | Denso IT Laboratory researcher's blog sites
こちらにまとめがあります。

ここでいう正規化定数とは、例えば指数型分布族の形を思い出すと
$\displaystyle p(x|\theta) = \frac{\exp(\eta(\theta) \cdot T(x))}{Z(\theta)}$
のような形をしていて、この分母の
$\displaystyle Z(\theta) := \int \exp(\eta(\theta) \cdot T(x)) dx$
を指します。（ベイズ統計の正規化定数とは違うもの）

単なる積分評価なのでモンテカルロ法やらの数値積分でも評価は可能ですが、計算時間や精度の問題があるため、積分によっては別のアプローチの方が良いことがあり、ホロノミック勾配法もそのうちの1つのアプローチ方法です。

ホロノミック勾配法では大きくは以下の2ステップで $\displaystyle Z(\theta)$ を評価します：

ステップ1. $\displaystyle Z(\theta)$ に関する微分方程式系を作る

ステップ2. 上記の微分方程式系を数値解析で解く

まずステップ1.の微分方程式系を作れる（存在する）かという疑問だが、 $\displaystyle Z(\theta)$ が名前の由来になっている"ホロノミック"であれば可能であることが知られています。

ではどんなものがホロノミックであるかというと

a. 有理函数

b. exp(有理函数), log(有理函数)

c. holonomic + holonomic, holonomic * holonomic

d. $\displaystyle \int \rm{holonomic}$

などがホロノミックであるので、統計で使われているほとんどのクラスはホロノミックと言えそう。

ここで、 $\displaystyle Z(\theta)$ がホロノミックであればグレブナー基底とかで機械的に微分方程式系を作ることが可能なことが知られています。たとえば、Macaulay 2にD-加群のアルゴリズムが入っているDmodulesパッケージがあるのでそれを使えば良いらしい。ただし、次元が高かったりすると計算量が大きくてかなりきついようで、手計算で微分方程式が出せるのであればそれを使うに越したことはない。

ステップ2.については初期値を適当に決めて（例えば級数展開をして気合で計算して）、Runge-Kutta 法などの数値解析手法で解くのですが、ここの部分の初期値や数値解析手法で時間や精度が決まってくると思うので、実務上はこの辺りが肝になってきそうである。初期値を決めた後はRのhgmパッケージでも計算できます。

では具体的な例を見てみよう。
方向統計学と呼ばれる方向にも興味がある時に使う分布で、von Mises-Fisher分布というものがあるのでそれを使ってみる。ここでは3変数von Mises-Fisher分布（on $\displaystyle S^2$ ）の分子だけ知っていたとしよう。つまり、
$\displaystyle p(x|\mu, \tau) = \frac{\exp(\tau \mu^T x)}{Z(\tau)}$
として $\displaystyle Z(\tau)$ をホロノミック勾配法で数値計算して求めてみよう。
ここで $\displaystyle \mu, \tau$ はそれぞれmean directionおよびconcentration parameterと呼ばれるもので、特に $\displaystyle \tau$ は以下の図のように $\displaystyle \tau =0$ の時は全ての方向に一様で、 $\displaystyle \tau = \infty$ の時はデルタ関数となる。

f:id:lambdahat:20180611154704p:plain — von Mises-Fisher分布のconcentration parameter

（図は http://people.csail.mit.edu/jstraub/download/straub2017vonMisesFisherInference.pdf より引用）

まずステップ1.については適当に微分や積分をすると、正規化定数Zに関して以下のような微分方程式を得ることができる：
$\displaystyle \frac{d^2}{d \tau^2}Z(\tau) = Z(\tau) - \frac{2}{\tau}\frac{d}{d \tau} Z(\tau)$

ステップ2.について、まず初期値を決めるためにベッセル関数を雑に近似させたものを使った。その後、上記の微分方程式と初期値をRのhgmパッケージのhgm.Rhgm関数に入れた。（hgmパッケージの計算例と実質同じ例です）

install.packages("hgm")
library(hgm)

#微分方程式を記述
dZ.fun = function(tau, Z, fn.params=NULL){
dZ = array(0, c(1, 2))
dZ[1,1] = Z[2]
dZ[1,2] = Z[1] - 2*Z[2]/tau
return(dZ)
}

#Zの初期値を決めるためにターゲットを雑に近似
Z.asymp = function(tau){
c(2*pi*exp(tau)/tau, 2*pi*(tau-1)*exp(tau)/tau^2)
}

#tauの適当な初期値とターゲットとしているtauを指定
tau_ini = 1
tau_target = 10

#Zの初期値を計算
Z_ini = Z.asymp(tau_ini)

#ホロノミック勾配法によるZ(10)の近似値
hgm.Rhgm(tau_ini, Z_ini, tau_target, dZ.fun)[1]

上記で計算された結果を見ると

[1] 13839.66

という数値が出てくる。
これがどれだけ正しいかを見るために、以下正確な値を計算しよう。

Z.exact = function(tau){
  c( 4*pi*(sinh(tau)/tau), 4*pi*(cosh(tau)/tau - sinh(tau)/tau^2) )
}

Z.exact(tau_target)[1]

この結果を見ると

[1] 13839.64

となっているので初期値が雑でも結構正確な値となっていることが分かる。

最後にこの方法が特異モデルにも使えるかが気になるところだが、微分方程式に特異点がある場合は何らかの対処が必要なようだ。ただ代数解析とも関係があるのでそのあたりもうまくやれそうな気もするが、かなりまじめに勉強をしないと分からなそうだということだけ分かった。

参考文献：
[1] JST CREST 日比チーム (2011), グレブナー道場, 共立出版.
[2] 日比, 竹村, 原, 東口, 清, グレブナー道場著者一同 (2015), グレブナー教室, 共立出版.

2018-05-17

純粋数学でアメリカ大学院受験をした時の昔話

もう10年前の話でTOEFLなどの試験をはじめ状況は全く違うと思うが、ネットで出てくる数学でアメリカ留学した先生達は優秀過ぎて参考にならないので、凡人の一記録を残しておくことで誰かの励みになってほしい。

大学や学部成績

・学部はニューヨーク州立大学バッファロー校という（少なくとも当時は）英語さえクリアしていれば誰でも入れる大学に1年から在籍していた（入ってみると日本人留学生はコミカレから編入する人達が大多数ということを知った）

・GPAは全体で3.74, 数学のみだと3.94 （1年目で取った歴史やライティングの成績がD(1.0)と壊滅的であった…。上述のような教養は簡単なコミカレで専門を4年制の大学で取る方がGPAを高く保てたであろう…）

・授業を取りまくって夏学期も取っていたので2年半で卒業できる状態になったが、大学院の入学のタイミングに合わせて卒業を半年遅らせて3年で卒業した

・よく知られているとおり日本人にとっては誰でも入れるアメリカの大学の学部数学の授業は凡人にとっても簡単なので、2年目からは大学院の授業を取っていた

・3年目の卒業研究はHeegaard Floer homologyの周辺をやっていた

・2年目の終わりか3年目の頭に数学科首席のawardを貰った

数学コンペ

・Putnam competitionというアメリカの大学生向けの数学オリンピックのようなものがあり、大学院を目指しているアメリカの数学科の人達は受けてる人が多かったと思う

・Top500に入ると各大学に名前のリストが送られるので、まずはtop500を目指す人が多い

・3年目の結果は296位/3545、1,2年目の結果は見つからないがそれぞれ550位と400位くらいだったと思う（3年目の結果は大学院の枠を争う同学年だけなら100位くらい？）

TOEFL

・当時はCBTという300点満点の試験でspeakingもなかった

・数学の大学院は250点以上を要求しているところがほとんどだった

・英語（というか興味が無い分野のお勉強）が苦手だったので253点というギリギリの点数だった

GRE

・数学の大学院だとGREはほとんど見られないという噂を信じほとんど対策をしなかった（することに越したことはない、特にsubject）

・当時はQuantitativeとVerbalは800点満点であった（当時と今での傾向の違いは知らない）

・Qantitative: 800, Verbal: 340, Writing: 4.0

・subjectの点数は覚えていないが、undergrad directorに点数悪すぎだろとイジられたのは覚えている（6時間で12問を解くPutnam competitionと違い、公務員試験のような短時間で大量の問題を解く試験で、このタイプの試験は今でも苦手である）

出願分野

・数学だと出願時点で純粋数学か応用数学かで枠が(大まかに)分かれていることが多い

・純粋数学の方が遥かにcompetitiveだが、純粋数学しか興味を持てなかったので純粋数学の方で出願した

・もちろん純粋数学がやりたくて純粋数学志望で入ったけど、入った後で応用数学の方に変える人も結構いる（逆は聞いたことがない）

出願結果(1年目)

・身の程知らずなのでTier 1のトップスクール（Harvard, MIT, Princeton, UC Berkeley, Stanford,...）のPh.D.プログラムだけに出願

・UPennだけwait listに入っていたが結局全滅した

その後…

・undergrad directorに全滅したので日本に帰って就職すると伝えたら、お前は数学を続けるべきだと説得され、授業料免除＋奨学金で通っていた大学のPh.D.プログラムに入れてもらった

・上記の説得の中で、この大学ではなくもっと良い大学院に行った方が良いと思うからTier 2くらいの再出願も来年やってみなと言われた

2年目の出願に向けて

・純粋数学だと1年目で論文書くなどしてアピールすることは不可能なので、入学早々Qualifying examを合格することを目標とした

・Ph.D.プログラム入学した月のQualifying examで無事合格した（学部の時から大学院の授業を取っていたため、この時点で申請すればMasterを取れる状態となった）

・勉強している中でK理論かその周辺の代数を研究したくなったので、出願先はそのあたりをやれそうなところだけにした

出願結果(2年目)

合格（全てPh.D.プログラムで学費免除+奨学金付+α）

・University of Southern California （ここに進学した）

・University of Utah

・Indiana University - Bloomington

・Northestern University

・University of Nebraska - Lincoln

Wait List

・UPenn （結局落ちた）

・Johns Hopkins University （USCの方が志望度が高く先に連絡が来ていたため、こちらから途中でお断りした）

余談

・上の学年と下の学年に同じく大学院を目指していた日本人の数学科の人がいたが、彼らも（純粋数学をやろうとしていたこともあり）大学院に全滅していたので、僕らのような平凡な人間だとアメリカの大学からでもアメリカのPh.D.プログラムに受かるのは簡単とは言えないと思う

最後に

・大変とはいえ、トップスクールでなくてもアメリカのPh.D.プログラムは刺激的で勉強になったので興味のある人は積極的に出願してみてほしい

・今調べてみたら数学でアメリカの大学院の留学について書かれた良い記事があったので、こちらの方が参考になるかと思います（私より遥かに優秀な人のケースではありますが…）

https://www.math.uci.edu/~takahasy/gradschool.html

2017-09-01

因子分析の情報量規準

とあるところで因子分析の情報量規準について見て、少し思うことがあったので簡単にまとめた。

因子分析のモデルは一般的には特異モデルであるため、BICの仮定が成立していない可能性がある。
（因子分析におけるBICの妥当性については後で述べるとして…）こういった場合、特異モデルの情報量規準として

WBIC (Watanabe, 2013)
sBIC (Drton & Plummer, 2017)

を使うことを検討することになるであろう。

因子分析に用いる場合、この2つでどちらの方が良さそうだろうか？
この2つの比較であれば以下の観点でsBICに軍配が上がると思われる。

対数周辺尤度の漸近展開の主要2項からのバイアス

因子分析であれば実対数閾値の理論値が知られている。したがってsBICの方が対数周辺尤度の漸近展開の主要2項はよく近似されると期待される。一方で、WBICはサンプルサイズが小さいときには対数周辺尤度の漸近展開の主要2項からのバイアスが小さくないことが知られている。
sBICに不利な点があるとするならば最適なパラメータのところをMLEで置き換えているため、それによるバイアスが入りうるところだと考えられる。

ばらつき

WBICはMCMCにより計算されるためばらつきが出てしまう。一方でsBICであればそのようなばらつきは出ない。

計算量

よく知られている通りMCMCの計算には非常に時間がかかってしまうためWBICの計算も大変になってしまう。sBICであればすぐに計算が可能である。

また、因子分析であればRのsBICパッケージがあるので簡単に計算できる。
実際に簡単なケースでシミュレーションをしてみよう。
シミュレーションの設定は、
Fitting a Bayesian Factor Analysis Model in Stan
を少し変えたものにしておく。具体的には無理やり特異的なケースにするため、
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3<- c(0.00, 0.00, 0.85, 0.80, 0.00, 0.75, 0.75, 0.00, 0.80, 0.80)
のl3をl2に似せて
l3 <- c(0.00, 0.95, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.35)
としておこう。こうすることで、実際は潜在因子が3つであるが2つに近い状態に見えてしまうであろう。
（因子分析的にこの設定に意味があるとは思えないが、元の設定からアドホックに特異なケースにしたものです）
簡単に設定をまとめておくと、真の因子は3つ（そのうちの2つは似てる）で、観測されている変数は10、サンプルサイズを50とし、因子数が1～5のうちどれが最も妥当であるかを周辺尤度を基準に決めようというもの。

WBICの計算は上のリンクのStanコードをちょっといじって計算させました（が、リンク元のコード部分が（ちょこちょこおかしなところがあり）正しい保証がないので、WBICの計算部分は今回は載せないことにしました）
後から気づきましたが、因子分析のstanコードは以下を参考にするのが良さそうです。
StanなMCMCで探索的因子分析やってみた | Kosugitti Labo ver.9

library(MASS)
library(sBIC)

D <- 3
P <- 10 
N <- 50

mu_theta <- rep(0,D) 
mu_epsilon <- rep(0,P) 
Phi <- diag(rep(1,D))
Psi <- diag(c(0.2079, 0.19, 0.1525, 0.20, 0.36, 0.1875, 0.1875, 1.00, 0.27, 0.27))
l1 <- c(0.99, 0.00, 0.25, 0.00, 0.80, 0.00, 0.50, 0.00, 0.00, 0.00)
l2 <- c(0.00, 0.90, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.30)
l3 <- c(0.00, 0.95, 0.25, 0.40, 0.00, 0.50, 0.00, 0.00, -0.30, -0.35)
L <- cbind(l1,l2,l3) 

fa <- FactorAnalyses(numCovariates = 10, maxNumFactors = 5)

iter <- 100

res <- data.frame(model_BIC=numeric(iter), model_sBIC=numeric(iter), model_WBIC=numeric(iter))

for(i in 1:iter){
  Theta <- mvrnorm(N, mu_theta, Phi) 
  Epsilon <- mvrnorm(N, mu_epsilon, Psi) 
  Y <- Theta%*%t(L) + Epsilon
  
  res.fa <- sBIC(Y, fa)
  res$model_sBIC[i] <- which.max(res.fa$sBIC)
  res$model_BIC[i] <- which.max(res.fa$BIC)

#WBIC部分は非掲載とした
}

summary(as.factor(res$model_sBIC))
summary(as.factor(res$model_BIC))

WBIC含めた結果は以下の通りとなった

手法＼推定された因子数	1	2	3 (true)	4	5
sBIC	0	0	80	19	1
WBIC	0	91	8	1	0
BIC	0	0	99	1	0

ということでWBICよりsBICの性能が高いことが分かる。
（WBICの結果の妥当性はやや不安が残るが、sBIC論文でのReduced Rank RegressionでのsBIC vs WBIC や他の因子分析のWBICの結果の報告とも傾向があっているので、傾向としては間違っていないと思われる）

しかし結果を見てみると、これら2つの特異モデル用の情報量規準よりBICの性能が非常に高い！

これについてはJRSS-BのsBIC論文のDiscussionでRafteryがLopes and West (2004)を引用して言及している通り "Good performance of standard BIC was also reported for factor analysis" なので、因子分析については実務的にはBICで良いのかもしれない。

sBICがそもそも何ものなのか、何をついて計算しているものなのかについてはまた今度書こうと思います。

参考文献
Drton M and Plummer M., A Bayesian information criterion for singular models, Journal of the Royal Statistical Society: Series B 2017, 79, 323--380

あざらしの日常

あざらしが腹を叩く片手間で書いたメモ