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

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

2016-05-21 2のd乗組み合わせデータ構造のための情報幾何

[][][][][][]Dually Flat Manifolds

  • 昨日、階層構造を持つ確率分布のための情報幾何について少し書いた
  • 特に、¥{0,1¥}^dなd元分割表のためのe-平坦、m-平坦なパラメタの取り方について、丁寧に確認してみたい
  • Rでパラメタ変換をしてみるのが、手っ取り早そうなので、それでやってみる
  • Rmd
---
title: "Dual flat parameterization for {0,1}^d"
author: "ryamada"
date: "Friday, May 20, 2016"
output: html_document
---

# Brief introduction

{0,1}の値を取るd個の変数$X_i;i=1,2,...$があるときに、その値セットの取り方($\mathbf{x} = (0,0,...,0),(1,0,...),...,(1,1,...,1)$)は
$D=2^d$通りある。

今、この$D$通りの取り方の生起確率ベクトルを$freq=(f_1,...,f_D)$とすれば、
$\sum f_i = 1$の制約があり、自由度は$D-1$である。


このような階層構造を持つ組み合わせについては、情報幾何的な観点から有用な、パラメタの取り方が2通り知られている。

It is know that there are tow useful ways to parameterize the frequency with the hierarchical structure.

この文書の目的は、頻度ベクトルを2通りのパラメタに置き換えたり、パラメタの値を指定して、頻度ベクトルを作ったりする道具を提供すること。

This document is written to offer utility functions to transform frequency vectors to two types of parameters and their inverse.




詳しい理論は ["Information Geometry on Hierarchy of Probability
Distributions" (IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 5, JULY 2001)](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwja3aaXwOXMAhXHXaYKHVSXCzAQFgghMAA&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F18%2F20133%2F00930911.pdf%3Farnumber%3D930911&usg=AFQjCNFd9avfcHWB3jmeVpfvnuWzCKo5rQ&sig2=OBOG6MXZgpPJvHd2JksQIQ) を参照のこと。

You can learn its theory in ["Information Geometry on Hierarchy of Probability
Distributions" (IEEE TRANSACTIONS ON INFORMATION THEORY, VOL. 47, NO. 5, JULY 2001)](https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwja3aaXwOXMAhXHXaYKHVSXCzAQFgghMAA&url=http%3A%2F%2Fieeexplore.ieee.org%2Fiel5%2F18%2F20133%2F00930911.pdf%3Farnumber%3D930911&usg=AFQjCNFd9avfcHWB3jmeVpfvnuWzCKo5rQ&sig2=OBOG6MXZgpPJvHd2JksQIQ)

## e-flat parameterization

$\mathbf{\theta}$は長さ$D-1$のパラメタベクトル。

$\mathbf{\theta}$ is a parameter vector with $D-1$ elements.

$\mathbf{\theta} = ((\theta_1,\theta_2,...,\theta_d,),(\theta_{1,2},...,\theta_{d-1,d}),...,(\theta_{1,2,...,d}))$

$\mathbf{x}=(x_1,...,x_d)$というパターンを取る確率
$$
\log{p(\mathbf{x})}=(x_1,...,x_d))= \sum \theta_i x_i + \sum \theta_{i,j} x_i x_j + ... + \theta_{1,2,...,d} x_1 x_2 ... x_d - \psi
$$

$psi$の値は$\mathbf{\theta}$が決まると(自由度がもうないので)一意に決まる。

The value of $\psi$ is not free because $\mathbf{\theta}$ has $D-1=df$ elements.

$\theta$ は1要素の寄与、要素ペアの寄与、要素トリオの寄与…に対応づく。


$\theta$ represents single effect, pairwise effect, triowise... and so on.

## m-flat parameterization

$\mathbf{\eta}$は、長さ$D-1$のパラメタベクトル。

$\mathbf{\eta}$ is a parameter vector with $D-1$ elements.

$\mathbf{\eta} = ((\eta_1,\eta_2,...,\eta_d,),(\eta_{1,2},...,\eta_{d-1,d}),...,(\eta_{1,2,...,d}))$

$\eta$の各要素は、変数(単独・組み合わせの)期待値を表す。(多元分割表で言えば)周辺頻度にも対応する。

$\eta$ corresponds to the expected value of single/combination of variables, or in other words, marginal frequency of multi-way table.

$$
\eta_{i} = Prob(x_i=1)\\
\eta_{i,j} = Prob(x_i=x_j=1)\\
...
\eta_{1,2,...,d} = Prob(x_1=x_2=...=x_d=1)
$$

# R codes in the ryamada22/Ronlyryamada github-package

# Utility functions and their usage example
```{r}
#library(devtools)
#install_github("ryamada22/Ronlyryamada")
library(Ronlyryamada)
example(my.haplotypes)
```

# Vignettes on functions

## Enumeration of $2^d$ patterns
d個の変数、$2^d$の場合の列挙。
1行目から順番に、値が1である変数の数が、0,1,...,dと増えるような並び方になっている。

$2^d$ patterns for d variables.
From the top row to the bottom, number of variables with value 1 increases like 0,1,...,d.


```{}
my.haplotypes <- function(d){
  n <- 2^d
	X <- matrix(0,d,n)
	cnt <- 1
	for(i in 1:d){
		tmp <- combn(d,i)
		for(j in 1:length(tmp[1,])){
			X[tmp[,j],cnt+j] <- 1
		}
		cnt <- cnt + length(tmp[1,])
	}
	return(X)
}
```
```{r}
d <- 3
my.haplotypes(d)
```

## Matrices to convert between frequency, theta and eta.

```{}
my.freq2eta.mat <- function(d){
  X <- my.haplotypes(d)
	n <- 2^d
	M <- matrix(0,n,n)

	for(i in 1:n){
		M[,i] <- apply(X >= X[,i],2,prod)
	}
	return(t(M[2:n,2:n]))
}
my.theta2freq.mat <- function(d){
  K <- my.freq2eta.mat(d)
	n <- 2^d
	M <- matrix(0,n,n)
	M[2:n,2:n] <- K
	M[,1] <- -1
	return(M)
}
my.freq2theta.mat <- function(d){
  M <- my.theta2freq.mat(d)
	return(solve(M))
}
# eta2freq.mat = inverse of freq2eta.mat
```
```{r}
my.freq2eta.mat(d)
my.theta2freq.mat(d)
my.freq2theta.mat(d)
```

## Transformations among three parameterizations

"freq", "$\theta$", "$eta$"3通りのパラメタ表現。

"$\theta$"の場合、自動的に決まるD番目の値を計算しておく方がよい。ほかの2つのパラメタシステムの場合は足して1になるところからの差、と言う形で簡単に計算できるが、$psi$はそうではなく面倒くさいから。


"freq" stands for regular frequency vector with $2^d=D$ elements.

"$theta$" stands for e-flat parameters.

"$eta$" stands for m-flat parameters.

For $\theta$, the residual parameter "\psi" should be also calculated because it is a bit difficult than other two systems' D-th value, even though it is automatically determined. 

```{}
my.theta2psi <- function(theta){
  len <- length(theta)
	d <- log(len+1,2)
	M <- my.theta2freq.mat(d)
	#tmp <- M %*% c(0,theta)
	tmp <- my.Inf.matrix.mult(M,c(0,theta))
	delta <- sum(exp(tmp))
	psi <- log(delta)
	return(list(psi=psi,M=M))
}
#' @export

my.theta2freq <- function(theta){
	out <- my.theta2psi(theta)
	psi.theta <- c(out$psi,theta)
	M <- out$M
	return(exp(my.Inf.matrix.mult(M,psi.theta)))
}
#' @export

my.freq2thetapsi <- function(freq){
	d <- log(length(freq),2)
	M <- my.freq2theta.mat(d)
	#return(M %*% log(freq))
	return(my.Inf.matrix.mult(M,log(freq)))
}
#' @export

my.freq2eta <- function(freq){
	d <- log(length(freq),2)
	K <- my.freq2eta.mat(d)
	return(K %*% freq[-1])
}
#' @export

my.eta2freq <- function(eta){
	d <- log(length(eta)+1,2)
	K.inv <- solve(my.freq2eta.mat(d))
	tmp <- K.inv %*% eta
	return(c(1-sum(tmp),tmp))
}
#' @export

my.theta2eta <- function(theta){
	freq <- my.theta2freq(theta)
	return(my.freq2eta(freq))
}
#' @export

my.eta2thetapsi <- function(eta){
	freq <- my.eta2freq(eta)
	return(my.freq2thetapsi(freq))
}
```

## Infinite values の扱い

freqベクトルに0が含まれると、$\theta$が無限大になる。

無限大を含めて計算できるようにするために、少し、工夫が必要。そのユーティリティ関数。

無限大と非無限大とを分離して行列計算し、最後にまとめる。

```{r}
my.matrix.Inf.separator <- function(M){
  if(!is.matrix(M)){
		M <- matrix(M,ncol=1)
	}
	M.inf <- M.noninf <- matrix(0,length(M[,1]),length(M[1,]))
	infs <- which(is.infinite(M))
	noninfs <- which(!is.infinite(M))
	M.noninf[noninfs] <- M[noninfs]
	M.inf[infs] <- sign(M[infs])
	return(list(M.noninf=M.noninf,M.inf=M.inf))
}

my.Inf.matrix.mult <- function(M1,M2){
	out1 <- my.matrix.Inf.separator(M1)
	out2 <- my.matrix.Inf.separator(M2)
	
	tmp.n.n <- out1$M.noninf %*% out2$M.noninf
	tmp.i.i <- out1$M.inf %*% out2$M.inf
	tmp.n.i <- out1$M.noninf %*% out2$M.inf
	tmp.i.n <- out1$M.inf %*% out2$M.noninf
	
	tmp.any.i <- sign(tmp.i.i + tmp.n.i + tmp.i.n)
	non.zero <- which(tmp.any.i!=0)
	tmp.any.i[non.zero] <- tmp.any.i[non.zero] * Inf
	tmp <- tmp.n.n + tmp.any.i
	return(tmp)
}
```

2016-05-20 Dually Flat Manifolds

[][][][][][]Dually Flat Manifolds

  • 資料
  • 確率分布を点としてもつ多様体を統計多様体という
  • 確率分布をn個のパラメタで表すことにすると、その統計多様体はn次元多様体
  • 多様体には多様体固有の特徴があって、それはリーマン計量だったりするわけで、局所座標系を入れる・入れないに関わらずその多様体の広がり具合というのは決まっている
  • 統計多様体では、座標系を入れる。パラメタで分布を表現することと同じ
  • パラメタを入れて、多様体全体をパラメタを使った座標表示にすると、格子が張り付いたことになる
  • この格子は、うまく張り付けると、ねじれたり捩れたりしないで素直に多様体上を流れていく
  • そんなよい感じの座標表示を"flat"と呼ぶ
    • 具体的には、3次元以上の多様体のときに、任意の3パラメタの間に延び広がり方・曲り方に偏りがないようなもののこと
    • 2通りのflatがあって、e-flat, m-flatと言うものがあり、この2つのflatは相互に密接な関係がある
    • またe-,m-flatsを表す式もよく似ている
    • e-flatの定義:確率分布が¥mathbf{¥theta}で表されている
      • E¥[¥frac{¥partial^2}{¥partial ¥theta_i ¥partial ¥theta_j}¥log{p(x,¥mathbf{¥theta})} ¥frac{¥partial}{¥partial ¥theta_k}¥log{p(x,¥mathbf{¥theta})}¥]=0
    • m-flatの定義:確率分布が¥mathbf{¥eta}で表されている
      • E¥[¥frac{1}{p(x,¥mathbf{¥eta})}¥frac{¥partial^2}{¥partial ¥eta_i ¥partial ¥eta_j}p(x,¥mathbf{¥theta}) ¥frac{¥partial}{¥partial ¥eta_k}¥log{p(x,¥mathbf{¥eta})}¥]=0
  • e-flatならm-flat、m-flatならe-flat
  • e-,f-両方のflatなパラメタ表現があるとよい(のだが、そうそう簡単に見つからないかもしれない)
  • ありがたいことに…
    • 多くの確率分布のパラメタ表現は、e-flatになるような表現で書きだせる。また、そのような表現を持つ確率分布を指数分布族と呼ぶ。そのような表現は
      • p(x,¥mathbf{¥theta})=e^{¥sum ¥theta_i k_i(x)- ¥psi (¥mathbf{¥theta})}
        • ただし、e-flatな分布関数が、すべて指数分布族なわけではない
    • また、m-flatになるような確率分布族も知られている。混合分布族
      • p(x,¥mathbf{¥eta})=¥sum ¥eta_i q_i(x) + (1-¥sum ¥eta_i) q_0(x)
        • ただし、m-flatな関数が、すべて混合分布族なわけでもない
  • e-flatなパラメタ表現がわかっても、m-flat表現がわからないとつまらないし、その逆も然り
    • ありがたいことに、ルジャンドル変換と呼ばれるルールで見つけられる。以下の変換ルールを見ればわかるように、¥mathbf{¥theta}¥mathbf{¥eta}とは、要素ごとに対応関係がある
      • ¥eta_i = ¥frac{¥partial}{¥partial ¥theta_i} ¥psi(¥mathbf{¥theta})= E¥[k_i(x)¥]
      • ¥theta_i = ¥frac{¥partial}{¥partial ¥eta_i} ¥phi(¥mathbf{¥eta})
        • ¥phi(¥mathbf{¥eta}) = E¥[¥log{p(x,¥mathbf{¥eta})}¥]
  • 対数尤度関数での関係
    • 確率分布を2通りのflatなパラメタ表現をすること、それらの相互関係について書いてきた
    • 確率分布関数を尤度関数とみなすとき、尤度関数の代わりに対数尤度関数を用いることも多い。確かに、上記の記載でもあちこちで対数化した関数が用いられていた
    • 以下では、対数尤度関数¥log{p(x,¥mathbf{¥theta})},¥log{p(x,¥mathbf{¥eta})}としてやり、そのパラメタ方向偏微分を成分とするベクトルを考えると、¥theta,¥etaとが「正規直交基底関係」のように見える(実際は、¥mathbf{¥theta}のそれぞれの成分同士は互いに直交ではなく、¥mathbf{¥eta}も同様なので、e-flat側とm-flat側との対応がきれいだ、ということになる
    • <e_i,e^*_j> = ¥delta_{ij}
  • Dual flatなパラメタ表記を用いると、次が言える
    • ¥psi(¥mathbf{¥theta}) + ¥phi(¥mathbf{¥eta}) + ¥mathbf{¥theta} ¥cdot ¥mathbf{¥eta} = 0
    • この式は、p(x,¥mathbf{¥theta}),p(x,¥mathbf{¥eta}という、「同じ関数」「多様体上の同じ点」の異なるパラメタ表記に関すること
    • ここで、「異なる点」についてこれを用いる
    • D¥[p:p’¥]=¥psi(¥mathbf{¥theta}) + ¥phi(¥mathbf{¥eta}’) + ¥mathbf{¥theta} ¥cdot ¥mathbf{¥eta}’というものが定義できて、これはp=p’のときに0でそうでないときには正の値をとる
    • これがKL divergenceとなっている
    • D¥[p:p’¥] = E_{¥mathbf{¥theta}}¥[¥log{¥frac{p(x,¥mathbf{¥theta})}{p(x,¥mathbf{¥theta}’)}}¥]
  • 関数間の遠近関係
    • 関数は統計多様体上の点なので、関数の違いを多様体上の「距離」として測りたい
    • しかしながら、KL divergenceの定義は、「自身」とのdivergenceを0となるような定義式を使って、「他者」のdivergenceを正の値にしており、2つの関数のどちらを基準にするかで値が変わる。そのことを受けて「距離」ではなく"divergence"と呼ぶ
    • e-,m-flatsという2通りのパラメタ系の良さは次の通り
      • e-系で表しているとき、e-flatなので、ごちゃごちゃ考えずとも、e-系のパラメタを使って比較的簡単にdivergenceを計量することができる
      • m-系の場合も同じ
      • 今、3つの関数p,p',p''があったときに、
      • p,p'について、e-系表現がわかっていれば、両者のdivergenceが簡単に測れる
      • さらに、p',p''について、m-系表現がわかっていれば、両者のdivergenceも簡単に測れる
      • p,p''の間のdivergenceは、片やe-系、片やm-系なので、簡単に測れないけれど、ここで、e-系とm-系とが「直交」していることが役に立ち、両者が単純に足し合わせればよいことが示せる
  • 関数間の遠近を統計量にする
    • KL divergence(対数尤度の差)の2倍の値はカイ二乗分布に従うことが知られており、尤度比検定に持ち込める
  • Flatなsubmanifolds
    • e-,m-flatsのよいところは、それぞれの座標系で、関数空間全体がベクトル空間のようになっていることである
    • 今、例えばe-系で表された関数空間全体のうち、e-系の線形和で表せる部分空間(submanifold)を取り出したとする
    • これはe-系表現があるからe-flatだし、e-flatならばm-flatではあるけれど、m-系で線形和に表せるかどうかというとそうはならない
      • e-系とm-系とは、相互に「直交」関係にあり、e-のすべての要素のそれぞれをm-系の全要素が分担して表現する、というような関係になっている(なっていがち)なために、e-の方の部分をとっても、m-の方は「個々の要素の寄与分」は小さくなっても、mの要素数を減らすわけにはいかないことなどを念頭に置くとわかりやすい
  • うまく辻褄の合うe-系submanifoldシリーズとm-系submanifoldシリーズ
    • e-系とm-系とのパラメタには、対応関係があるので、それについて注意をしつつ、パラメタを除いていくと(パラメタの値を固定していくと)、入れ子になったe系submanifoldsシリーズとm-系submanifoldsシリーズができる。それは、直交性を利用して活用することができる

2013-05-09 ぱらぱらめくる『情報幾何の基礎概念』

[][]メモ

  • 多くの見慣れた・聞き慣れた確率分布は指数分布族のメンバーである
  • こちらがWikiで、こちらが日本語のサイト
  • どうして同じ族に属し、また「exponential」と冠されるかと言えば、このメンバーの確率密度関数がf_X(x|¥theta) = h(x)exp(¥eta (¥theta) T(x)-A(¥eta))と表されるから
  • Wikiの表に沿って、これをやってみよう
    • 2項分布の確率密度分布は¥begin{pmatrix} n¥¥x¥end{pmatrix}p^x(1-p)^{n-x}
    • f_X(x|¥theta) = h(x)exp(¥eta (¥theta) T(x)-A(¥theta)):こちらを使うと、Wikiの表から、
    • h(x)=¥begin{pmatrix}n¥¥x¥end{pmatrix}: Base measure
    • ¥eta(¥theta) = ¥log{¥frac{p}{1-p}}: Natural parameter(s)
    • T(x) = x: Sufficient statistic
    • A(¥eta) = -n¥log{1-p}: Log-partition A(¥theta)
    • f_X(x|¥theta) = ¥begin{pmatrix}n¥¥x¥end{pmatrix}exp( ¥log{¥frac{p}{1-p}}x+n¥log{1-p})となるが、これは
    • f_X(x|¥theta) =¥begin{pmatrix}n¥¥x¥end{pmatrix}(¥frac{p}{1-p})^x(1-p)^nとなって、確かに同じ
  • nカテゴリあって、その確率の和が1になるようなものも、式変形することで、指数分布族に含まれることがわかる。ちらの第3ページ