2018-03-17

確率質量分布で表すRandom Exchangeable Partitions

Kingmanのpaintbox〜単位線分タイリングがrandom exchangeable partitionsの表現であることがわかったが、実際、どんなrandom exchangeale partitionsが現れるのかは確率事象なので
確率事象として生成されるときにどんな確率事象なのか、そのときにどのような分布になるのか、ということが気にかかる
前の記事では有限個線分の個数を指定し、一様ディリクレ乱数によって線分の長さを定めて生成した。この方法だと、有限個線分の指定個数がかなり大きな制約条件とはなっているが、これも、一つの生成確率事象ルールではある
以下では、いくつかの方法を記述する
どんな条件があるかと言えば、 $\sum_{i=0}^\infty p_i=1;p_i \ge 0$ 。確率が0より大とすれば確実に無限個に分割。確率が0のものが無限個あるとすれば、有限個への分割もこの定義でカバーされる
自然数・非負整数を大とする確率質量分布は使える
- 有限個の項に対応する多項分布の確率質量に比例した単位線分分割は、Kingmanのタイリングの無限小分割を持たないバージョンとなる

k <- 10
sort(rdirichlet(1,rep(1,k)))

- ポアソン分布の確率質量は0,1,...,無限大のそれぞれに正の確率質量が与えられその和は1

plot(dpois(0:100,3),type="h")
plot(sort(dpois(0:100,3)),type="h")

- ポアソン分布よりも分散が大きい確率質量分布に負の二項分布がある

plot(dnbinom(0:100,3,0.3),type="h")
plot(sort(dnbinom(0:100,3,0.3)),type="h")

- ディリクレ分布、ポアソン分布、負の二項分布のいずれも、ガンマ分布を使って次のように表現できる
  - k項のディリクレ乱数は、尺度パラメターが1のガンマ分布に従うk個の乱数 $\Gamma_1,...,\Gamma_k$ を発生し $\frac{\Gamma_i}{\sum_{i=1}^k \Gamma_i}$ とすることで生成できる。このとき $\sum_{i=1}^k \Gamma_i$ は有限の値を取ることが知られている。このときのガンマ分布の形状パラメタベクトルがディリクレ分布の形状パラメタになっている(らしい)
  - ポアソン分布の場合は1個のガンマ分布乱数をその期待値として指定することでガンマ分布と紐づく
  - 負の二項分布の場合は、ガンマ乱数を発生させ、その乱数に対応するポアソン分布から乱数を一つ取り出すことを繰り返すとその結果が負の二項分布となると知られています。負の二項分布のsizeとprobとのパラメタとガンマ分布のshape,rateパラメタとの関係はshape = size, rate = prob/(1-prob)です)こちら)

N <- 10^4
a <- c(0.1,0.2,0.3)
Rdirichlet <- rdirichlet(N,a)
Rgamma <- cbind(rgamma(N,a[1],1),rgamma(N,a[2],1),rgamma(N,a[3],1))
Rgamma <- Rgamma/apply(Rgamma,1,sum)
plot(sort(Rdirichlet[,1]),sort(Rgamma[,1]))

N <- 10^5
sz <- 4
pr <- 0.3
sp <- sz
rt <- pr/(1-pr)
Rnb <- rnbinom(N,sz,pr)
RgammaPois <- rep(0,N)
for(i in 1:N){
	tmp <- rgamma(1,sp,rt)
	RgammaPois[i] <- rpois(1,tmp)
}

plot(sort(Rnb),sort(RgammaPois))

2018-03-17

Kingman's theorem、Random Exchangeable Partitions

Random Exchangeable Partitions Poisson Dirichlet過程 Coalescent Size-biased Pick Kingman

無限大()のRandom Exchangeable Partitions
- - 限大にするとちょっと厄介
こんな方法(KingmanのPaintboxの方法)というのがある
- 1,2,...,Nという数列を長さ１の線分に見立てて、それを分割する
- ただし、Nは無限大なのでこの線分上には無限個の自然数が並んでいる
- これを有限個に区切って、そこには線分の長さに比例した数の自然数が帰属するとする。ただし、全体が無限個の自然数に対応するから、有限の長さの線分はどんなに短くても無限個の自然数が対応する
- このようにすると、無限個の自然数を有限個にしか分割できない
分割数を無限にする
- 無限個にある分割が有限の長さに対応するとそれを合算すると長さ１の線分に納まらない
- なので、無限個の分割には無限小の長さを対応付ける
- これを実現するために、長さ１の線分を有限個に分割し、そのうちの１つは、自然数が１個帰属する小さな分割の集まりとする。それ以外の分割線分には長さに応じた無限個の自然数が帰属するとする
- こうすることで、無限個の $\frac{1}{\infty}$ の長さの線分の集合があることを表現できる
- Exchangeableであることを考慮すれば、長さ１の線分を有限個に適当に分割し、無限分割用の線分を取り除ける。それ以外の線分は長さの順にならべても一般性を失わない
- そんな風にして無限個への分割を作成し
- さらに、無限個の自然数を割り当てることにする。ただし、無限個の自然数の割り当ては無限に続けなくてはならないので、適当な数までランダムに割り付けることにする

# Kingman's paint box

k <- 10
library(MCMCpack)
r <- rdirichlet(1,rep(1,k))
r. <- cumsum(r)

n <- 10^5
u <- runif(n)
my.label <- function(x){
	tmp <- which(r. > x)
	tmp[1]	
}
lb <- unlist(sapply(u,my.label))

plot(r,c(table(lb)))
abline(0,n,col=2)

このKingmanのPaintboxの方法は単位線分のタイリングとも称されるが、exchangeable random partitionとKingmanのpaintbox タイリングとは1対1対応することが知られている

2018-03-15

ポアソン・ディリクレ過程とハプロタイプ頻度

ハプロタイプポアソン・ディリクレ過程ポアソン過程ディリクレ分布

サンプリングして有限個数の多項分布を生み出すことに関するものと、長さ１を分割する・無限分割することに関するものに大きく２分され、その２つが相互に関連する、という構図になっている
中華料理店過程とその一般化としてのピットマン・ヨー過程
多項標本
分割
- Pitman’s model of random partitions
- ピットマン確率分割と関連する話題

2018-03-06

Priorを調整しながらベイズ

ベイズ事前分布

あるpriorでMCMC ベイズを回して事後分布を得るとする
そのpriorは、ある基準で選ばれたpriorだが、別の基準だと「変数変換」しないといけないとする
そんなpriorの重み変換をすることができるのか、できるならどうやるのかの調べもの
こちらは、この用に適しているか不明ながら、メモ

2018-02-28

FDR: Benjamini-Hockberg

FDR Knockoff Benjamini-Hochberg

昨日の記事はKnockoff 変数を用いたFDRの制御の話
FDRといえば、Benjamini-Hochbergもある
これは、「ある閾値で変数の取捨選択をするとする」ときに、すべの変数が帰無仮説ＯＫ変数だったとした場合に、何個の変数がFalselyに帰無仮説を棄却するかの期待個数を算出する
その数を、実際に帰無仮説を棄却した変数の数で割ってやる。この中には、真に関連のある変数が入っているので、全部帰無の仮定の下での期待個数より大きいと考えられる
この割合をある値q(たとえば0.05)以下にしてくれるような閾値のうち、もっとも小さい値をとろう、という戦略
そうすることで、その割合(FDR)をq以下にコントロールしつつ、閾値を越える変数の数を最大にできる、という作戦
p個の変数が独立であることを前提に考えている
以下、Rでやってみる
- p個の変数があり、k個は真に関連があり、残りは帰無
- 統計量が正規乱数で得られるとする(ほかの統計量でももちろんよい)。p値でやっても、もちろんよい
- 閾値をどうするかを探索するためのプロットが黒と赤水平線
- 実際に偽陽性の割合を緑でプロット
- 黒カーブが赤水平線と交わるところを閾値にすることになり、そこより右側では、緑は赤水平線より下に来る。これがFDRをコントロールしてある、という意味

p <- 10000
k <- 200

z <- rnorm(p)
for(i in 1:k){
	z[i] <- rnorm(1,rnorm(5,2),1)*sample(c(-1,1),1)
}

t <- seq(from=0,to=5,length=1000)

q <- 0.05

x <- rep(0,length(t))
y <- rep(0,length(t)) # no. selected features
w <- rep(0,length(t)) # no. true positive
for(i in 1:length(x)){
	x[i] <- p * pnorm(t[i],lower.tail=FALSE)*2/max(sum(abs(z)>=t[i]),1)
	
	y[i] <- sum(abs(z)>=t[i])
	tmp <- which(abs(z)>=t[i])
	w[i] <- y[i]-sum(tmp %in% 1:k)
	
}

plot(t,x)
abline(h=q,col=2)
y. <- y
y.[which(y==0)] <- 1
points(t,w/y.,col=3)

2018-02-27

Knockoff 変数

Knockoff

Knockoff 変数を使ったFDRについての概説記事はこちら
説明変数 Xがnxp行列(nサンプル、p個の説明変数)であるときに
$\begin{pmatrix}X,\tilde{X}\end{pmatrix}^T\begin{pmatrix}X,\tilde{X}\end{pmatrix} =\begin{pmatrix}\Sigma, \Sigma - diag(s)\\ \Sigma-diag(s), \Sigma \end{pmatrix}$
この形は2p個の変数の分散共分散行列になっており、Positive definite
したがって、 $diag(s)$ (対角行列)の取り方に制約がある
その制約がある中で、なるべく、 $\Sigma - diag(s)$ の対角成分を0に近づけたい。なぜなら、オリジナル変数と対応するKnockoff 変数との関係はできるだけ「無関係〜対角成分が0」であってほしいから
ではどんな制約があるか、というと、Schur complementを使って以下のようにする
- $G= \begin{pmatrix}A,B\\C,D \end{pmatrix}$ であるとき、 $G/D = A - B D^{-1}C$ である。今、 $A= D = \Sigma,B = C=\Sigma- diag(s)$ であるから
- $G/\Sigma = \Sigma -(\Sigma - diag(s)) \Sigma^{-1} (\Sigma - diag(s))$ となり、式変形していくと
- $G=diag(s) (2 \Sigma - diag(s))$ となる。両辺がpositive definiteなので、右辺の２項 $diag(s)$ と $\Sigma -diag(s)$ との両方がpositive definite
この条件を満たす $\diag(s)$ を探せばよい
Knockoff パッケージでは、対角行列 $diag(s)$ の対角成分を全部同じにして探す方法(create_equicollerated)を採用していたり、凸最適化であるSemidefinite optimization(こちらやこちらにて探索する方法を採用していたりする
このような方法でKnockoff 変数が作れるのは、 $n \ge 2p$ の場合
それ以外の場合は、 $\begin{pmatrix}X,\tilde{X}\end{pmatrix}$ の条件を満足するように工夫しながら、 $\tilde{X}$ を算出する

2018-02-27

Knockoff 変数によるFDR

Knockoff FDR PPV

資料はこれ(基本)とこれ(GWAS等への拡張)
Rのパッケージはknockoffで、そのgithubがこちら
考え方の基本
- FDRをしたい
  - 多変量解析をしていて、いくつかの変量は従属変量に意味のある寄与があり、残りの変量は意味がない、というように振り分けたい
  - その振り分けにあたって、意味があるとみなした変量の中に含まれる偽陽性の割合(FDR)をコントロールしたい
- Knockoff 変数を使う
  - p個の説明変数があるときに、ある方法でp個のKnockoff 変数を設定する。このKnockoff変数は、p個の説明変数のそれぞれに対応する変数として定める
  - Knockoff 変数は次のような変数である
    - Knockoff 変数p個の相互関係は、オリジナルのp個の変数の相互関係と「同じ」である。同じ、というのは、分散共分散行列が同じという意味で「同じ相互関係」である
    - あるKnockoff変数と、対応するオリジナル変数以外のオリジナル変数との関係は、「対応するオリジナル変数とそれ以外のオリジナル変数との関係」と同じである
    - あるKnockoff変数と、対応するオリジナル変数との関係は、できるだけ「無関係」である
  - このようなKnockoff 変数をどのようにして作成するか、はさておき(その説明はこちらに別掲)、このようなKnockoff 変数があったとき、従属変数と無関係な説明変数とその対応Knockoff 変数とは、どちらも従属変数と無関係であり、従属変数と関係のある説明変数の場合、オリジナルの方は従属変数と関係するが、Knockoff 変数の方は、オリジナル変数と「できるだけ無関係」に取られているわけであり、従属変数とは無関係である
- Knockoff 変数を用いたFDRのコントロール
  - オリジナル変数とKnockoff 変数とを合わせて2p個の説明変数で多変量解析をすることを考える
  - 無関係なオリジナル変数とそのKnockoff 変数とでは、どちらが「採択的に優先される」かは、等確率と期待される
  - 関係のあるオリジナル変数とそのKnockoff 変数とでは、「オリジナルの方」が「採択的に優先される」はずである
  - この関係を利用すると、2p個の変数に「採択ルール」を適用し、ある基準を越えたら採択することにする。ただし、オリジナル変数がその対応するKnockoff 変数より「優先」されている場合に限るとすると、結局、p個のオリジナル変数の(添え字の)どれかを採択していることになる
  - ただし、このようにして採択された(オリジナル)変数の中には、「たまたま基準を超え」、かつ、「たまたまオリジナルの方が対応するKnockoff 変数より優先されていた」ものが含まれている
  - どの(オリジナル)変数が、「たまたま採択された」のかは、不明だが、そのようにして採択される変数の個数は、だいたい合っている。なぜなら、変数の数をpから2pへと倍にして、その上で、「たまたまオリジナルが上という確率0.5」で選択してあるから
  - ただし、この「たまたま」の個数はこれだけではわからない。しかしながら、Knockoff 変数が「採択ルール基準を越え」、それがオリジナル変数よりも優先度が高くなっているようなものがどれかはわかるし、その数は数えられる。このようなものの中には、稀に、「本当にオリジナル変数が意味のある変数である」こともあるかもしれないが、基本的には、このようになるのは、従属変数と関係のないオリジナル変数とそのKnockoff 変数のペアに限るし、「たまたまKnockoff の方がオリジナルよりも優先される確率は0.5」なので、この個数がFalse discoveryの個数と同じとみなせる
  - したがって、「採択基準を超え、かつ、Knockoff 優先になっている個数」を「採択基準を超え、かつ、オリジナル優先になっている個数」で割った値が、FDR
どんな「採択ルール・優先度」を使うか
- うまく使えるルール・優先度〜統計量にはある制約があるが、それはおいておいて
- たとえば、LASSOをやったときに、ペナルティの重み変数をどこまで小さくするとモデルに組み込まれるか、という値は、ここで言う「優先度」として使える
- また、Stepwise selectionをするときの、組み込まれる順番も優先度として使える