Hierarchical Models とは「階層型モデル」と訳せる。階層型モデルは例えば
あるコイン工場にある偏りを ω 、そこで作られた複数のコインの偏りを θs
コインの偏りは、工場に特有の偏りに依存する(その逆はない)。つまり、工場の偏りは階層の上位で、下位のコインの偏りに影響するという関係。これが(単純化した)階層型モデル。
θs が ω に依存していることを表したのが次の式 9.1
p(θ, ω | D) ∝ p(D | θ, ω) p(θ, ω)
= p(D | θ) p(θ | ω) p(ω) ... (9.1)
左辺は Bayes' rule 、右辺は Bayes' rule の分子で、両辺の比例関係を示している。ここでは推定するパラメータが θ, ω 。
この式の二行目で一行目の右辺を変換する。尤度 p(D | θ, ω) は、データは θ に依存するので ω はデータと無関係なので、p(D | θ) に置換される。事前確率 p(θ, ω) は、p(θ | ω) = p(θ, ω) / p(ω) の関係から p(θ | ω) p(ω) と置換される。
この変換結果から、データは θ に依存( p(D | θ) )し、θ は ω に依存( p(θ | ω) )することがわかる。先のコイン工場の例では
The estimate of the bias of any one coin depends on the estimate of factory bias, which is influenced by the data from all the coins.
工場の偏りの推定値 ω に、コインの偏り θs が依存している。そして、工場の偏りは全コインのデータを反映する。
注意すべきは「観測データ(全コインのデータ)に影響されるのは、θ ではなく ω」で、「データは θs に依存する」が「その逆ではない」ということ。ランダムウォークでサンプリングするのは θs だが、その θs が依存するのは ω で、そしてデータが影響するのは ω という関係。
文章だけはわかりにくかもしれない。この階層型モデルを学びながら連想したのは 因子分析 。「数学と理科が得意の背景には、理数系という共通の因子がある」(個人的には「理数系」「文系」のカテゴライズは嫌いだが...)、などの因子である。θ というサンプリングデータの結果の背後には ω という因子、つまり「工場特有の偏り」があるということ。
この辺の考察は今後の課題とするが、非常に興味深い。
A Single Coin From a Single Mint
このタイトル、「1つの造幣局から1種類のコイン」という意味。語呂が良いので、以降は英文のままで使う。
左は図 9.1 で、本書で用いられる階層型モデルの表記方法。パラメータ関係の分かりやすさ、JAGS コードとの容易な対応など、本表記のメリットは多い。
この図が表すモデルは、例えば、ある種類のコインの偏りを θ 、そのサンプリングデータを yi 、θ が依存する上位の変数が ω 、これで "A Single Coin From a Single Ming"。
以下、詳細は省くが、各階層を数式で表す。
ω ~ dbeta(a, b) ... (9.3)
a = ω(k - 2) + 1, b = (1 - ω)(k - 2)+1
ω = (a - 1) / (a + b -2)
k = a + b
θ ~ dbeta(ω(k-2)+1, (1-ω)(k-2)+1) ... (9.4)
注意:式 9.4 において、図 9.1 では定数 K だが、現実的なモデルでは変数 k 。詳細は以降の投稿で扱う。
p(ω) = beta(ω | Aω, Bω) ... (9.5)
ω, k と θ の関係
θ の ω への依存は式 9.4 から分かるが、concentration k との関係は
the magnitude of k is an expression of our prior certainty regarding the dependence of θ on ω
k は、θ の ω への依存性について、事前の確実性を表現したもの
より具体的な表現では
The estimated credible values of ω will be similar to θ if K is large, but will be much less certain if K is small.
推定した ω の値は、 K が大きければ θ に近似し、K が小さければ不確定となる。
以下の R の実験においても ω = 0.25, k = 1000, 100, 10, 5 の場合、k が大きほど θ の平均値は 0.25 に近似している。
> theta <- function(k,omega) rbeta(100,omega*(k-2)+1,(1-omega)*(k-2)+1)
> for (i in c(1000,100,10,5)) print(summary((theta(i,0.25))))
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.2148 0.2466 0.2527 0.2522 0.2596 0.2845
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1224 0.2207 0.2586 0.2543 0.2876 0.3634
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.03385 0.17530 0.29610 0.30210 0.40050 0.62990
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.05589 0.23180 0.36660 0.38090 0.49340 0.88880
「視覚化」に続く。

0 件のコメント:
コメントを投稿