ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2016-09-20 インポータンス・サンプリング、Exponential tilting 指数型関数族

[][][]インポータンス・サンプリング、Exponential tilting 指数型関数族

  • n人の投擲選手が居て、その人たちの投擲記録について賭けをすることを考える
  • n人の投擲選手の記録はガンマ分布に従うと予想する。ガンマ分布を選ぶのは、0以上の1ピークの分布であって、まあまあいい感じのモデルだから
  • ある選手の記録は以下のような確率分布になる
x <- seq(from=0,to=100,length=101)
m1 <- 30
sc1 <- 5
sh1 <- m1/sc1

y1 <- dgamma(x,shape=sh1,scale=sc1)
plot(x,y1,type="l")

f:id:ryamada22:20160920122903j:image

  • 選手人数n=4として、それらの誰が投げるかは等確率で決まるものとする
  • このとき、投擲記録の確率密度分布は以下のようになる
n <- 4
m <- c(20,40,50,60)
sc <- c(1,2,1,3)
sh <- m/sc

ys <- matrix(0,n,length(x))
for(i in 1:n){
	ys[i,] <- dgamma(x,shape=sh[i],scale=sc[i])
}
par(mfcol=c(1,2))
matplot(x,t(ys),type="l")
y <- apply(ys,2,mean)
plot(x,y,type="l")
par(mfcol=c(1,1))

f:id:ryamada22:20160920124236j:image

  • さて、賭けというのは、こういう賭けである
  • 賭ける人は、平均Mで標準偏差Sの正規分布を指定して良いことになっている。4選手の誰かが投げて出た記録に対して、指定した正規分布の確率密度を得点として得られる、というのが、この賭けのルールである
  • 今から賭けようと思うが、正規分布のパラメタの組(M,S)について、得点の期待値が知りたい。さて、どうするか
  • 地道にやってみよう
  • xの値を0から小刻みに増やして、それぞれの生起確率を出して、そのときの得点を積み上げていくことにする
    • それぞれのチェックポイントには生起確率と得点がある
M <- 40
S <- 20

n.r <- 10^4 # チェックポイント数
x.check <- seq(from=0,to=1000,length=n.r)
#x.check <- seq(from=0,to=1000,by=0.01)

Point <- 0
sum.prob <- 0
for(i in 1:length(x.check)){
	tmp.prob <- 0
	for(j in 1:n){
		tmp.prob <- tmp.prob + dgamma(x.check[i],shape=sh[j],scale=sc[j])
	}
	tmp.prob <- tmp.prob/n
	Point <- Point + dnorm(x.check[i],M,S) * tmp.prob
	sum.prob <- sum.prob + tmp.prob
}
Point <- Point/sum.prob
Point
> Point
[1] 0.01479605
    • これは、言い方を変えると、ある指定範囲に一様乱数のようなもの(等間隔指定点)を発生させて、その生起確率で重みをつけて、重み付き平均を出している
    • 重み付き平均を期待値とするのはインポータンスサンプリングのやり方だから、「一様分布」と「確率密度分布」とで重みを決めて、乱点のようなものとして、等間隔指定点セットを取っている場合に相当する
  • 生起確率密度分布から乱数を発生させて推定してみる
rs <- matrix(0,n,n.r/n)
for(i in 1:n){
	rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
}

hist(rs)

Point <- 0
for(i in 1:length(rs)){
	Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
}

Point
> Point
[1] 0.01485484
    • これは、まったく重みを付けずに、超地道に乱数を発生させて期待値を出している
  • さて。インポータンスサンプリングでは、先に一様乱数のようなものを使ったが、もっと別の確率分布を使ってみよう、というもの
    • どうしてそうするとよいか、というと、重み付き平均を出すときの、「重み」x「値」が大きめな値で揃っているときに、期待値の推定精度が上がることが知られているから
    • なので、何かうまい分布を取りたいな、という話
    • 例えば、投擲距離の理論分布の平均と標準偏差を、上で生成した乱数のそれを取って、その平均と標準偏差を持つような正規乱数を発生させ、その正規分布とガンマ分布の混合分布である理論分布との比を重みづけにしてやってみる

mm <- mean(rs)
ss <- sd(rs)

rs <- rnorm(n.r,mm,ss)

Point <- 0
for(i in 1:length(rs)){
	tmp.prob <- 0
	for(j in 1:n){
		tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
	}
	tmp.prob <- tmp.prob/n
	tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
	tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
	Point <- Point + tmp.point * tmp.weight
}
Point
> Point
[1] 0.01480977
  • ここまでで、変な確率密度分布をする確率変数があって、それに応じて、何かしらの得点を想定し、その期待値を求める話をしてきた
    • 区間からのの等間隔点を出すのもやったが、その代わりに一様乱数を使うことにして、同様のことをやってみる
    • ただし、乱点の数を少な目にして、モンテカルロ推定値がばらつくようにして、そのバラツキが大きいか小さいかを見てみる
M <- 40
S <- 20

n.r <- 100

n.iter <- 1000

Points.unif <- Points.gamma <- Points.norm <- rep(0,n.iter)
for(ii in 1:n.iter){
	#x.check <- seq(from=0,to=1000,length=n.r)
	#x.check <- seq(from=0,to=1000,by=0.01)
	rs <- runif(n.r,0,1000)
	Point <- 0
	sum.prob <- 0
	for(i in 1:length(rs)){
		tmp.prob <- 0
		for(j in 1:n){
			tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
		}
		tmp.prob <- tmp.prob/n
		tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
		tmp.weight <- tmp.prob/dunif(rs[i],0,1000)
		Point <- Point + tmp.point * tmp.weight
	}

	Points.unif[ii] <- Point
	rs <- matrix(0,n,n.r/n)
	for(i in 1:n){
		rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
	}

	Point <- 0
	for(i in 1:length(rs)){
		Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
	}

	Points.gamma[ii] <- Point

	rs <- rnorm(n.r,mm,ss)

	Point <- 0
	for(i in 1:length(rs)){
		tmp.prob <- 0
		for(j in 1:n){
			tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
		}
		tmp.prob <- tmp.prob/n
		tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
		tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
		Point <- Point + tmp.point * tmp.weight
	}
	Points.norm[ii] <- Point

}

Points.three <- cbind(Points.unif,Points.gamma,Points.norm)

boxplot(Points.three)

> apply(Points.three,2,mean)
 Points.unif Points.gamma  Points.norm 
  0.01498932   0.01479045   0.01477595 
> apply(Points.three,2,sd)
 Points.unif Points.gamma  Points.norm 
0.0068368258 0.0003667666 0.0006546649 
    • わかることは、生起確率密度そのものから乱数が発生できる場合(中央)はばらつきが小さく、一様乱数の場合はばらつきが大きいこと、正規乱数は結構、ばらつきが小さいことである
    • このことから、一様乱数より正規乱数の方がよいことである。それは、正規乱数が一様乱数よりも、本当の生起確率分布と得点値との積の分布に近いから

f:id:ryamada22:20160920135527j:image

  • では、生起確率密度が高いときに得点が小さいような場合はどうなるだろうか
    • このときに正規乱数の方を、「賭ける得点分布」にしてやってみる。なぜなら、正規乱数の密度x得点に近いときが良いのだから、『得点が高い場合には、投擲記録の出る確率に大差がないような場合』には、得点の分布のみの高低が影響しそうだから
M <- 150
S <- 30
mm <- M
ss <- S
n.iter <- 1000

Points.unif <- Points.gamma <- Points.norm <- rep(0,n.iter)
for(ii in 1:n.iter){
	#x.check <- seq(from=0,to=1000,length=n.r)
	#x.check <- seq(from=0,to=1000,by=0.01)
	rs <- runif(n.r,0,1000)
	Point <- 0
	sum.prob <- 0
	for(i in 1:length(rs)){
		tmp.prob <- 0
		for(j in 1:n){
			tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
		}
		tmp.prob <- tmp.prob/n
		tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
		tmp.weight <- tmp.prob/dunif(rs[i],0,1000)
		Point <- Point + tmp.point * tmp.weight
	}

	Points.unif[ii] <- Point
	rs <- matrix(0,n,n.r/n)
	for(i in 1:n){
		rs[i,] <- rgamma(n.r/n,shape=sh[i],scale=sc[i])
	}

	Point <- 0
	for(i in 1:length(rs)){
		Point <- Point + dnorm(rs[i],M,S) * 1/length(rs)
	}

	Points.gamma[ii] <- Point

	rs <- rnorm(n.r,mm,ss)

	Point <- 0
	for(i in 1:length(rs)){
		tmp.prob <- 0
		for(j in 1:n){
			tmp.prob <- tmp.prob + dgamma(rs[i],shape=sh[j],scale=sc[j])
		}
		tmp.prob <- tmp.prob/n
		tmp.point <- dnorm(rs[i],M,S) * 1/length(rs)
		tmp.weight <- tmp.prob/dnorm(rs[i],mm,ss)
		Point <- Point + tmp.point * tmp.weight
	}
	Points.norm[ii] <- Point

}

Points.three <- cbind(Points.unif,Points.gamma,Points.norm)

boxplot(Points.three)

> apply(Points.three,2,mean)
 Points.unif Points.gamma  Points.norm 
6.670910e-06 6.483460e-06 6.389868e-06 
> apply(Points.three,2,sd)
 Points.unif Points.gamma  Points.norm 
2.819655e-06 9.335974e-06 7.721868e-06
    • 今度は、正規乱数の方がばらつきがわずかだが小さくなったことが解る
  • このようなことを勘案して、どんな乱数を使うのがよいかを検討してサンプリングするのがインポータンス・サンプリング
  • ただし、生起確率密度がわかっていて、そのあたりの得点が高いときには、正規確率密度から乱数を発生させられるものなら発生させればよい。しかしながら、生起確率密度はわかるけれど、そこからの乱数発生関数がないような場合には、近くて乱数発生関数があるような乱数で代用し、その代り、重みづけについては、密度の比を取るのがよい
  • また、生起確率密度からの乱数発生ができる場合であっても、大きな得点を得られるところの密度が低いと、乱数が無駄になる。そんなときには大きな得点が得られるところに高密度な別の乱数発生分布を設定して、それとの比を重みにするのがよい。インポータンスサンプリングの主の使い道はこれだろう(か?)
  • それ以外の場合で、「得点を表した乱数」を発生することが良いような場合であっても、重みを付けるには、本当の生起確率密度が計算できなければ比が出せないことは同じ
  • さて。Exponential tilting
    • これは、生起確率分布が指数型分布族に属していて、乱数発生に適切な分布も指数型分布族に属しているときには、両者を指数型分布族表現してやり、その比(それがインポータンス・サンプリングの重み)が指数表現になっていて簡単だよ、という話
  • 参考:こちらとか、こちらとか、こちらとか

2016-09-17 あらためて指数型分布族・情報幾何

[][]あらためて指数型分布族・情報幾何

  • 指数型分布族のことを整理するためにはいろいろな段階があると思う
  • 現時点では、(全部で何段階かあるのかわからないけれど、そのうちの)第3段階あたりに居るような気がする
  • この先、階段を上るためにも、自分の各段階での理解を書いておく
  • 第3段階だけを書いておくと、(しばらく離れてしまって、自分で読み返したときに)第3段階までも到達できないかもしれないという危惧があるから

[][]私の第0段階〜あらためて指数型分布族・情報幾何

  • まだ指数型分布族という名前を知らない
  • 指数分布という名前は知っているかもしれない
  • 確率分布・尤度に関しては、以下のようなことを知っている
  • 確率分布はパラメタを使った関数。積分すると1
  • 関数が分布の「様子」を決め、パラメタがその縮尺などを決める
  • 確率密度関数があると、ある値を観測する確率が計算できる
  • 逆に、ある値が観察されているときに、その確率分布に照らして、パラメタの値に関して尤度が計算できる
  • 同じ式が確率も表し、尤度も表す
  • 尤度関数という呼び方もある
  • 尤度を最大にするパラメタの値を最尤推定値と言う
  • 最尤推定値を求めるには、尤度関数の微分が0のところを探せばよい。尤度関数の対数を取って、対数尤度関数の微分が0のところを探すのもよい
  • 指数関数と対数関数がある。掛け算と足し算との関係にある。確率事象が相互に独立に起きるときは確率を掛け合わせればよい。掛け算が面倒くさければ確率の対数を足し合わせればよい

[][]私の第1段階〜あらためて指数型分布族・情報幾何

  • 指数型分布族という名前に出会う
  • 以下のようなことを理解する
  • 非常に多くの有名な確率分布がすべて指数型分布族に属する
    • たまに指数型分布族に属さない分布もある。二つの正規分布の混合や、コーシー分布とか。混合正規分布は普通の統計解析アプローチで面倒くさいやつだ、統計的学習やBUGS/STAN学習とかを使う必要があった。コーシー分布は分散が定義されていない、とか何がどう変なのかわからないけれど、異端な分布だった。このような特徴と「指数型分布族」に属さないことに関連があるのか〜とぼんやり思う
  • 多数の分布(離散も連続も、Xがスカラーでもベクトルでも)が指数関数を使った統一された式で表せる
    • その例はこちらにびっしりと書いてある
    • 分布関数はパラメタを用いて表現されていたが、そのパラメタを変えることで指数関数の中に押し込むのが指数型分布族のやりかた
  • 分布関数というのは、「分布の本体」と「その形を実際に決めるパラメタ」とからなる。したがって分布関数はパラメタが作る空間に連続的に広がっている。その広がっている空間ではパラメタの取り方を変えることができる
  • 指数型分布族の定義式を見る
    • パラメタのみの項、確率変数のみの項、パラメタと確率変数とが関係しあっている項の3つの項からなっている
  • 指数型分布族の性質を見る
    • そのうちの1つは、自分が指数型分布族に興味を持ったきっかけなので、意義くらいは分かる
    • それ以外にたくさんの性質があって、どれも、「どうしてそうなるのか」を考えるとくらくらする
    • そもそも、色々な性質の説明文の中に、聞いたことはあるけれど、よくわかっていない単語が立て続けに出てきて、途方に暮れる
  • この段階での理解をまとめると
  • 『たくさんの大事な分布関数は、パラメタの取り方を変えることで、指数関数表現ができる。そのことを通じて、分布関数というものを統一的に理解することができて、とても有意義らしいが、この世界は記載方法があっちこっちでまちまちだし、素敵な諸性質も難しい』

[][]私の第2段階〜あらためて指数型分布族・情報幾何

  • 指数型分布族の性質や諸要素について、色々な呼称・断片的な知識に曝露されつつも、それらの有機的な関係が乏しい段階
  • 指数型分布族の定義
    • 定義は1つだが、式の書き方はそうではない
      • すべての項を指数の肩に乗せる書き方と、パラメタ単独項を前に出す書き方の2法が主流である
      • パラメタと確率変数とが関係しあっている項では、¥theta ¥cdot T(x)のように、パラメタはそのまま、確率変数は関数となっているが、場合によっては、パラメタの方も、関数にする(U(¥theta) ¥cdot T(x))ことがある
  • 指数型分布族の性質はごちゃごちゃしている
    • 共役事前分布を考えやすい
      • お互いに指数型分布族のメンバーなので、「本質的には同じもの同士」
    • 最尤推定しやすい
      • 尤度関数が指数型なので対数尤度関数にすればただの足し算。そこにパラメタと変数との線形関係があるだけ
    • 一般化線形回帰が説明できる
      • 指数の肩に線形式が入っていることとかんけいがある?
    • 正規化項は積率母関数、キュムラント生成関数
      • 積率母関数は確率変数を指数の肩にのせた関数の期待値を関数だが、この「指数の肩」が「指数型」とうまく折り合う
    • 正規化項から凸とわかる
      • 尤度を考えパラメタの最尤推定をするときには、正規化項の微分が大事。正規化項は確率密度の積分を1にするための項なので、正
    • 情報幾何で使う
      • 観測データが与えられるとパラメタが作る関数の多様体が見える
    • フィッシャー情報量
      • 情報多様体の局所曲率のこと。最尤探索しやすいところは曲率が高いところ
    • 接続と双対平坦
      • 多様体は局所(の座標の取り方)の他に、隣り合う局所同士の座標の接続具合を定めることで一意になる。その取り方に双対関係にあるよいものがある
    • 十分統計量
      • 指数表現にすると、確率変数の「あるべき姿」が見える
    • 自然パラメタ
      • パラメタの取り方はいろいろにできるが、指数表現という「すばらしい表現」のときのパラメタにはそれなりの意味がある

[][]私の第3段階〜あらためて指数型分布族・情報幾何

  • ごちゃごちゃしてきて収集がつかない
  • これをうまく整理するには、やはり数式が必要だし、納得するには式変形を見てみる必要がある
  • 数式とその変形を追いかけると、個々の性質の「意味」がどんどんかすんでいくので、そのあたりの折り合いをどうつけるかが課題
  • 現在、この段階の真っただ中なので、うまく書けるか不安だが、書いてみる
  • 確率密度関数・確率マス関数、パラメタ
    • 確率変数¥mathbf{x}
    • パラメタは取り方によっていろいろある
      • ¥mathbf{¥theta}=(¥theta_i),¥mathbf{¥eta},¥mathbf{¥xi}の3つとする。ただし、そのパラメタにも共通するように書きたいときには、「一般的なパラメタ」という意味で¥alphaを使うことにする
        • ¥thetaは指数型表現でのパラメタ
        • ¥etaは情報幾何で言うところの¥thetaの双対にあたるパラメタ
        • ¥xiは指数型にしないで、個々の分布関数の特徴をよく表したパラメタとする
        • ¥theta,¥eta,¥xiはそれぞれ行き合える。¥xi(¥theta)と書いたら、それはパラメタ¥xi¥thetaの関数として表した、ということ
  • 指数型関数の定義
    • ¥mathbf{¥theta}を使う
    • p(¥mathbf{x}|¥mathbf{¥theta}) = exp(¥mathbf{¥theta} ¥cdot T(¥mathbf{x}) - A(¥mathbf{¥theta}) - B(¥mathbf{x}))
      • これを少し違った形にすることもある
      • p(¥mathbf{x}|¥mathbf{¥theta}) = h(¥mathbf{x})g(¥mathbf{¥theta}) exp(¥mathbf{¥theta} ¥cdot T(¥mathbf{x}))
        • h(¥mathbf{x}) = exp(-B(¥mathbf{x})),g(¥mathbf{¥theta}) = exp(-A(¥mathbf{¥theta}))
        • log(h(¥mathbf{x})) = -B(¥mathbf{x}),log(g(¥mathbf{¥theta})) = -A(¥mathbf{¥theta})
  • 確率密度(マス)関数として扱う
    • パラメタが与えられたものとして、指数型関数を¥mathbf{x}の関数と見る
      • ¥int_{¥mathbf{x}} p(¥mathbf{x}|¥mathbf{¥theta}) dx = 1は確率密度関数の定義
      • exp(A(¥mathbf{¥theta})) = ¥int_{¥mathbf{x}} exp(¥mathbf{¥theta} ¥cdot T(¥mathbf{x}) - B(¥mathbf{x}))となるから¥frac{1}{g(¥mathbf{¥theta})}=exp(-A(¥mathbf{¥theta}))は正規化項(で正の値)
      • 確率変数の値の様子について知ることができる。様子と言ったらモーメント(積率)。積率母関数は、指数型表現の¥mathbf{x}によらない成分が確定する。積率母関数とその対数版であるキュムラント母関数は、それぞれ非¥mathbf{x}依存成分の関数を指数として見るか、指数関数の外に出すかに対応する
        • g(¥mathbf{¥theta}) = exp(-A(¥mathbf{¥theta}))は積率母関数、log(g(¥mathbf{¥theta})) = -A(¥mathbf{¥theta})はキュムラント母関数
        • この積率母関数はT_i(¥mathbf{x})の積率を教えてくれる。したがってT_1(x)=xのような単変量のときには、いわゆる確率変数そのもののモーメントがわかることになる
        • こちらの記事にて確かめられる
        • 確率変数があったときに、パラメタのみの変更で分布の位置や縮尺のようなものが変わるが、そのときにモーメントも変わるはず。であれば、そのようなモーメントに関する情報は、パラメタのみが支配する項によらなければならないのではないか。そんなことを考えると、A(¥mathbf{¥theta})がモーメントに関する情報を担っていることは納得が行きやすい。そのうえで、パラメタに関する偏微分とその次数がモーメントに逐一対応する、というのは、¥mathbf{¥theta}が特別なパラメタであるがための恩恵に相当すると思っておけばよさそう
      • 指数型表現は、ある意味で特別な意味を持つ表現であることがわかった。¥mathbf{x}のみの項が積率母関数(キュムラント母関数)を表しているのがその例。では、T(¥mathbf{x})は?これは、¥mathbf{x}の取り方を、特別なパラメタである¥mathbf{¥theta}に合わせて設定した、¥mathbf{x}の『あるべき姿』。¥mathbf{¥theta} ¥cdot T(¥mathbf{x})は内積だが、内積というのは線形汎関数と見ることもできるように、特別な関係であって、双対的な関係にある(こちら)
      • このことから、指数型分布ではないコーシー分布の場合、積率母関数が(指数型分布のようには少なくとも)書けない、ということがわかる
  • 尤度関数として使う
    • パラメタが与えられたものとして、指数型関数を¥mathbf{x}の関数と見る
      • 対数尤度関数はlog(p(¥mathbf{x}|¥mathbf{¥theta})) = ¥mathbf{¥theta} ¥cdot T(¥mathbf{x}) - A(¥mathbf{¥theta}) - B(¥mathbf{x})
      • 最尤推定は微分して0なので¥frac{¥partial}{¥partial ¥theta_i}log(p(¥mathbf{x}|¥mathbf{¥theta})) = 0
      • したがってT_i(¥mathbf{x}) - ¥frac{ ¥partial A(¥mathbf{¥theta})}{¥partial ¥theta_i}=0を解くのだが、これは簡単
        • このことから、指数型に書けない混合正規分布のような統計モデルのときには、簡単に最尤推定ができなくて、何か別のことをする必要が出ることもわかる
      • 尤度関数として見る、とは、¥mathbf{x}は定数扱いにするということ。すると対数尤度関数は¥mathbf{¥theta} ¥cdot T - A(¥mathbf{¥theta}) - C¥mathbf{¥theta}が張る空間での多次元曲面になっており、これが対数尤度関数の多様体
      • 多様体には、局所の曲率と局所に張り付ける座標の連なり具合(接続)を考えるのが、多様体的発想。局所曲率をフィッシャー情報量(情報行列)と言い、それは多様体にとってのリーマン曲率行列
      • 共役分布。指数型分布を1個以上の独立な観察について積み重ねることは、尤度関数を掛け合わせること。指数型関数の掛け合わせは相変わらず指数型。同じ関数を尤度関数と見るか確率密度関数と見るかで、パラメタを推定するための関数と見ることもできれば確率変数を発生させるための関数と見ることもできる。この味方の変換関係が共役な関係にある2種類の指数型分布
      • 二重平坦座標系。情報幾何では指数型関数の対数尤度関数が作る情報多様体では¥mathbf{¥theta}がe-平坦座標系をなしていることが知られている。情報多様体の上にうまく測地線を乗せた接続のことである。他方、e-平坦座標系と双対関係にあるm-平坦座標系というものがある。これは¥mathbf{¥eta}で表され、¥eta_i = E(T_i(¥mathbf{x}))の関係にある
      • 情報幾何ではエントロピーを気にする。エントロピーはlog(p)の平均値なので¥int p log(p)のこと。いずれにせよ、log(p)が大事なわけであるが、指数型分布族では、この対数確率・対数尤度が¥mathbf{¥theta} ¥cdot T(¥mathbf{x}) - A(¥mathbf{¥theta}) - B(¥mathbf{x})という単純な形をしている。このように表現できることから、log(p)の多様体が¥thetaに関してアフィンになっている(こともわかるという)。実際には多数の観察がもたらす尤度を問題にするのだが、log(p)ベースで行うときは、ただの加算(アフィン)になる
      • 情報幾何では、指数型分布族にも混合型分布族にも、離散にも連続にもあてはまるルールがある(この記事では、指数型分布族の表現が情報幾何でのe平坦座標を考えるのに有用だったので、両方を扱ってきた)
        • n個の値を取りうる確率分布を考える。連続分布の場合にはnは無限大。そのときある値を基準にして、そのほかの値の生起確率を基準値の生起確率に対する比でとり、その対数をとったものが¥mathbf{¥theta}パラメタ、生起確率そのものを取ったものを¥mathbf{¥eta}パラメタと取る。エントロピーが定義できて、KLダイバージェンスもとれる。KLダイバージェンスに対してピタゴラスの定理が成立し、射影も定義できる。エントロピー関数の2階微分がフィッシャー情報量になる(こちら)
      • ただし、指数型分布族はやはり好ましい性質を持っており、指数関数型で表した時のパラメタ¥mathbf{¥theta}を使うと、対数尤度関数が平坦になっている(e-平坦)、また、いつも双対平坦になっており¥eta_i = E(T_i(¥mathbf{x}))がそれ
      • 検定は情報多様体上の観測点(それはフルモデルの中にある)を部分多様体(部分モデルの設定の仕方によって変わってくる)上の点に対応付けること。それは、一次近似をするだけなら、局所の接空間の話になり、その接空間は尤度関数全体が決めるから、部分モデルの取り方によらず、精度は一定。もし、それよりも細かいことを問題にすると曲面を考えることになるが、曲面を考えるためには曲率が必要で、曲率を考えるときには接続を考慮しないと、観測点と部分も出る点との遠近関係が言えなくなるし、接続の取り方によって、最近点は変わってくる。同じ観測データに対して、実際、2次の近似をするときには、その誤差の遠近は、e曲率とm曲率と、部分も出るのパラメタの取り方が産む曲がり具合(曲率)とに分解できる。e曲率は共通だし、m曲率も知りたいことは同じなので、共通で、最後のパラメタの取り方による曲率というのが、検定手法依存の項になる。ちなみに最尤推定というのは、この部分モデルのパラメタの取り方による曲率をm曲率に一致させる方法なのだという。世の中で通用している検定は、1次の近似までは同じで、ばらつきが一緒。2次の項は方法によって、どういうときに大きくなりがちで、どういうときに小さくなりがちかに差が出てくる。それは手法の得手不得手の情報幾何学的説明となる

2016-09-14 t検定の一般化としてのホテリングのt-squared

[][][]t検定の一般化としてのホテリングのt-squared

  • 1次元正規分布を「原点からの距離」で考え直すと、-xとxとが同じ距離にあることになり、0 <= xに関する分布を考えることとなり、カイ分布になる。距離の二乗で考えることにすればカイ二乗分布
  • d次元正規分布を「原点からの距離」で考え直すと、原点から等距離にある点を集約して考えることになり、非負の距離値に関して考えれば、自由度dのカイ分布になるし、距離の二乗で考えれば自由度dのカイ二乗分布になる
  • t分布は、1次元に正規分布で分布している2群の平均値の差の検定などに使われる分布。これは1次元の世界のもの
  • 多次元空間に正規分布があるとする。2群あって、平均が同じか違うかに興味があるときには、多次元正規分布の平均値(平均座標)の異同を検定したい
  • そんなときに「多次元平均座標の差」をスカラー値として測定して、検定に持ち込みたい。t検定を多次元版に一般化したい、と言い換えてもよい
  • それがホテリングのt-squared(二乗しているので非負値を取る)
  • 一般化になっていることを納得するには、一般次元で動いているホテリングのt-squaredの1次元版とt検定とが対応していることを見ればよい
  • やってみる
library(MCMCpack) # 多次元正規分布を作るのに逆Wishart分布にて行列を作るため
library(Hotelling) # ホテリング検定のため
df <- 3 # 次元
N1 <- 1000 # 第1群サンプル数
N2 <- 1200 # 第2群サンプル数

n.iter <- 1000 # 帰無仮説の下でデータを発生させる回数
H <- p <- rep(0,n.iter)
for(i in 1:n.iter){
	S <- riwish(df,diag(rep(1,df)))
	X1 <- rmvnorm(N1,mean=rep(0,df),sigma=S)
	X2 <- rmvnorm(N2,mean=rep(0,df),sigma=S)
	tmp <- hotelling.test(X1,X2)
	H[i] <- tmp[[1]][1]
	p[i] <- tmp[[2]]
}
# 統計量とp値を見てみる。p値は一様分布になっている
par(mfcol=c(1,2))
hist(unlist(H))
hist(p)
par(mfcol=c(1,1))

# 次元を1にして、t検定の結果と比べる
df <- 1
N1 <- 1000
N2 <- 1200

n.iter <- 1000
H <- p <- t <- tp <- rep(0,n.iter)
for(i in 1:n.iter){
	S <- riwish(df,diag(rep(1,df)))
	X1 <- rmvnorm(N1,mean=rep(0,df),sigma=S)
	X2 <- rmvnorm(N2,mean=rep(0,df),sigma=S)
	tmp <- hotelling.test(X1,X2)
	H[i] <- tmp[[1]][1]
	p[i] <- tmp[[2]]
	tmp2 <- t.test(X1,X2)
	t[i] <- tmp2[[1]]
	tp[i] <- tmp2[[3]]
}
par(mfcol=c(1,2))
# t値とH's t-sqの平方根は絶対値が同じ
plot(sqrt(unlist(H)),t)
plot(p,tp)
par(mfcol=c(1,1))

2016-09-13 Stanで多次元正規分布

[][][]Stanで多次元正規分布

  • 多次元データが取れる場合と、多次元だけれど、その一部次元のデータしか取れない場合っていうのは、どういう違いがあるのか…
  • multinormal2.stan
data {
  int N_subjects;
  int N_items;
  matrix[N_subjects,N_items] y;
}

parameters {
  vector[N_items] mu;
  cov_matrix[N_items] Sigma;

}

transformed parameters {
}

model {
  matrix[N_items,N_items] identity;
  mu ~ normal(0,100);
  identity <- diag_matrix(rep_vector(1.0,N_items));
  Sigma ~ inv_wishart(N_items,identity);
  for(i in 1:N_subjects)
    y[i] ~ multi_normal(mu,Sigma);
}
  • multinormal4.stan
data {
  int N_subjects[3];
  int N_items;
  matrix[N_subjects[1],2] y12;
  matrix[N_subjects[2],2] y23;
  matrix[N_subjects[3],2] y31;
}

parameters {
  vector[N_items] mu;
  cov_matrix[N_items] Sigma;

}

transformed parameters {
  cov_matrix[2] Sigma12;
  cov_matrix[2] Sigma23;
  cov_matrix[2] Sigma31;
  vector[2] mu12;
  vector[2] mu23;
  vector[2] mu31;
  
  Sigma12[1,1] <- Sigma[1,1];
  Sigma12[1,2] <- Sigma[1,2];
  Sigma12[2,1] <- Sigma[2,1];
  Sigma12[2,2] <- Sigma[2,2];
  
  Sigma23[1,1] <- Sigma[2,2];
  Sigma23[1,2] <- Sigma[2,3];
  Sigma23[2,1] <- Sigma[3,2];
  Sigma23[2,2] <- Sigma[3,3];
  
  Sigma31[1,1] <- Sigma[3,3];
  Sigma31[1,2] <- Sigma[3,1];
  Sigma31[2,1] <- Sigma[1,3];
  Sigma31[2,2] <- Sigma[1,1];
  
  mu12[1] <- mu[1];
  mu12[2] <- mu[2];
  mu23[1] <- mu[2];
  mu23[2] <- mu[3];
  mu31[1] <- mu[3];
  mu31[2] <- mu[1];
}

model {
  matrix[N_items,N_items] identity;

  mu ~ normal(0,100);
  identity <- diag_matrix(rep_vector(1.0,N_items));
  Sigma ~ inv_wishart(N_items,identity);
  
  for(i in 1:N_subjects[1]){
    y12[i] ~ multi_normal(mu12,Sigma12);
  }
  for(i in 1:N_subjects[2]){
    y23[i] ~ multi_normal(mu23,Sigma23);
  }
  for(i in 1:N_subjects[3]){
    y31[i] ~ multi_normal(mu31,Sigma31);
  }
}
library(mvtnorm)
library(rstan)
library(MCMCpack)
N_subjects <- rep(100,3)
N_items <- 3

S <- riwish(N_items,diag(rep(1,N_items)))

mu <- rnorm(N_items)
tmp <- rmvnorm(N_subjects[1],mean=mu,sigma=S)
y12 <- tmp[,1:2]
tmp <- rmvnorm(N_subjects[2],mean=mu,sigma=S)
y23 <- tmp[,2:3]
tmp <- rmvnorm(N_subjects[3],mean=mu,sigma=S)
y31 <- tmp[,c(3,1)]

dat1 <- list(N_subjects=N_subjects,N_items=N_items,y12=y12,y23=y23,y31=y31)

out1 <- stan(file="multinormal4.stan",data=dat1)

ymp <- rmvnorm(sum(N_subjects)/3*3,mean=mu,sigma=S)
dat2 <- list(N_subjects=N_subjects[3],N_items=N_items,y=tmp)
out2 <- stan(file="multinormal2.stan",data=dat2)

2016-09-11 勝手にStan

[][]勝手にStan

  • Stanを使えば事後分布が乱択的に得られる
  • けれど、単純なモデルではつまらないし、普通に回帰等をする方がよっぽどよい
  • どれくらいモデルが複雑だとStanがよいかを言い切るほどよくわかっていないけれど、階層ベイズはそうなのだろう(こちら)
  • SNPの遺伝因子評価に使うとどうなるか、いじってみる
  • ひとつは、性・年齢が単純でなく関わっているとき。ただしStanでのモデル化は単純なので、結局、glm()と大差ないはず
  • もうひとつは、SNP-phenotype関連と、phenotype 発現量関連と、SNP-eqtl関連を統合するようなもので、phenotypeが入っている関連では個人が共通で、eqtlは全く別の個人セットとする(とおもったけれど、もう少し複雑にしないとやりたいことをやっていないようだった…要修正)
  • ひとつめ
functions{}

data{
	int<lower=0> nsample;
	int<lower=0,upper=1> phenotype[nsample];
	int<lower=0,upper=2> homo[nsample];
	int<lower=0,upper=1> hetero[nsample];
	int<lower=0,upper=150> age[nsample];
	int<lower=0,upper=1> sex[nsample];
}

transformed data{
}

parameters{
	real rhomo;
	real<lower=0,upper=1> rhetero;
	real rsex;
	real rage;
}

transformed parameters{
	vector[nsample] risk;
	for(i in 1:nsample){
		risk[i] <- rhomo * homo[i] + rhetero * rhomo * hetero[i] + rsex * sex[i] + rage * age[i];
	}
}

model{
	rhomo ~ normal(0,1000);
	rhetero ~ uniform(0,1);
	rsex ~ normal(0,1000);
	rage ~ normal(0,1000);
	
	for(i in 1:nsample){
		phenotype[i] ~ bernoulli(inv_logit(risk[i]));
	}
}

generated quantities{}

gp <- c(0.49,0.42,0.09)
n <- 10^3
g <- sample(0:2,n,replace=TRUE,prob=gp)
homo <- hetero <- rep(0,n)
homo[which(g==2)] <- 1
hetero[which(g==1)] <- 1
sex <- sample(0:1,n,replace=TRUE,prob=c(0.3,0.7))
library(MCMCpack)
age <- sample(20:70,n,replace=TRUE,prob=rdirichlet(1,rep(1,length(20:70))))
age. <- age-40
age.[which(age.<0)] <- 0

ph <- (g+1)^2 + sex * 2 + (age./max(age.))^1.8 * 5 + ((age-30)/max(age-30))^2*(-3.4)
ph. <- ph + rnorm(n,0,mean(ph))

ph.obs <- (ph.-min(ph.))/(max(ph.-min(ph.)))
ph.obs. <- ph.obs > mean(ph.obs)
glm(ph.obs. ~ g,family="binomial")

PH <- as.numeric(ph.obs.)

dat <- list(phenotype=PH,nsample=n,age=age,sex=sex,homo=homo,hetero=hetero)
stan.out <- stan(file="SNPcovariate.stan",data=dat)

stan_hist(stan.out,pars=c("rhomo","rhetero","rage","rsex"))
glm.out <- glm(PH~g+hetero+age+sex,family="binomial")

  • ふたつめ
functions{}

data{
	int<lower=0> nsample;
	int<lower=0,upper=1> phenotype[nsample];
	int<lower=0,upper=2> homo[nsample];
	int<lower=0,upper=1> hetero[nsample];
	int<lower=0,upper=150> age[nsample];
	int<lower=0,upper=1> sex[nsample];
}

transformed data{
}

parameters{
	real rhomo;
	real<lower=0,upper=1> rhetero;
	real rsex;
	real rage;
}

transformed parameters{
	vector[nsample] risk;
	for(i in 1:nsample){
		risk[i] <- rhomo * homo[i] + rhetero * rhomo * hetero[i] + rsex * sex[i] + rage * age[i];
	}
}

model{
	rhomo ~ normal(0,1000);
	rhetero ~ uniform(0,1);
	rsex ~ normal(0,1000);
	rage ~ normal(0,1000);
	
	for(i in 1:nsample){
		phenotype[i] ~ bernoulli(inv_logit(risk[i]));
	}
}

generated quantities{}

ngenotype <- 1000
g <- sample(0:2,ngenotype,replace=TRUE,prob=c(0.49,42,0.09))
pre.ph <- g * 3 + 1 + rnorm(ngenotype,1)
ph <- as.numeric(pre.ph > mean(pre.ph))

nexpression <- 500
gexpression <- sample(0:2,nexpression,replace=TRUE)
expression <- rep(0,nexpression)
for(i in 1:nexpression){
	expression[i] <- rnorm(1,gexpression[i],1)
}

blood <- rep(0,ngenotype)
for(i in 1:ngenotype){
	blood[i] <- rnorm(1,g[i],1)
}

dat <- list(ngenotype=ngenotype,neqtl=nexpression,genotype=g,phenotype=ph,expression=expression,gexpression=gexpression,blool=blood)

stan.out <- stan(file="GEBP.stan",data=dat)