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

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

2017-03-18 HLA領域メモ

[][]HLA領域メモ

  • Tevfik DORAK’s homepage
  • その彼のよいサマリー
  • まとめ方。上記を参考にして
    • (どうしてそんなに複雑なのか?)
    • HLA領域(extended MHC)の位置、大きさ
    • 関連Phenotypes:免疫系、非免疫系
    • HLAの複雑さの特徴は何か?
      • 遺伝子の数、密度
      • その内訳:HLA,non-HLA
      • extended LD
    • HLA領域の免疫系と非免疫系
    • タイプ分け(DNA側からと免疫認識側から)、歴史的・進化的側面
    • 複雑さとフェノタイプ関連の検出機構(どうして劣性評価するか、Extended haplotypeの存在と、組み合わせ効果の検出効率)

2017-03-12 へいきん

ryamada222017-03-12

[][][]へいきん

  • こちら,こちら,こちら,こちらでやっていることのこと
  • 観測値の標本平均値は¥frac{¥sum_{i=1}^n x_i}{n}
  • これはf(x|m) = ¥sum_{i=1}^n (x_i-m)^2を最小にする¥hat{m}
    • ¥frac{d f}{dm} = ¥sum_{i=1}^n 2(x_i-m) = 2 (¥sum_{i=1}^n x_i -n ¥times m) = 0の解
    • 実数直線にx^2という距離の定義を入れ、そのWasserstein distanceの和を最小にするようなm ( ¥hat{m} = argmin_m ¥sum_{i=1}^n W(m,x_i)^2_2)
  • 別の見方をすると、平均m(分散はなんでもOK)なる正規分布からの標本であるとみなしたときの尤度・対数尤度を最大にするm
  • 「へいきん」というのが、算術平均(足して割る)なので、対数尤度(尤度は足し合わすことができる)を「距離コスト」としてWasserstein距離の下での「へいきん」が最尤推定値であり、標本平均
  • こう考えれば、正規分布ではない統計モデルの下でのモデルパラメタの最尤推定は、パラメタが張る空間に尤度で距離コストを定めたうえでの、Wasserstein 平均を求めることが、モデルパラメタ最尤推定となる
  • 逆に、標本が複数あって、それらが指し示す「ベストな推定」をする、とは、空間と距離の定義を入れた上でWasserstei distanceの総和の最小化問題になる
  • 推定課題は、空間の置き方とそこの距離の置き方のバリエーション
  • さて…。
  • また、http://d.hatena.ne.jp/ryamada/20170310では、この最小化推定をするアルゴリズムとして、空間を全探索する代わりにエントロピーを使って1次元多様体に沿った探索に落としている
  • その落とし方は、正則化項(LASSOとか、圧縮とか)とよく似た式に変わることに注意
  • こちらこちらもメモ

2017-03-09 統計遺伝学分野で生き抜くための数学〜精選19

ryamada222017-03-09

[][][]KLダイバージェンス

  • 双対平坦空間では、拡張ピタゴラスの定理が成り立つので、斜辺を共有する直角三角形が、2つの円周(球面)を構成する

f:id:ryamada22:20170310071652p:image

---
title: "分布の違いを測る Measure dissimilarity between distributions KLダイバージェンス KL-divergence"
author: "ryamada"
date: "2017年3月6日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 平坦と直線と直交 Flatness, straight lines and perpendicular

## 直線 Straigt lines

- $\theta,\eta$座標系はそれぞれ平坦座標系であるから、$\theta_i$一定、$\eta_j$一定の点(上図の4種類の線)は、「直線」である。Because both $\theta$ and $\eta$ coordinate systems are flat, points with same value of $\theta_i$ or $\eta_j$ are in "straight lines", that are shown in the figure above with 4 colors.

- $\theta_1,\theta_2$の線形和で表される線も「直線」。$\eta_1,\eta_2$の線形和で表される線も「直線」。Points that are expressed as linear sum of $\theta_1$ and $\theta_2$ are in "straight lines" as well. Same for $\eta_1$ and $eta_2$.

## 直交

- $\theta$座標系と$\eta$座標系は双直交であるから、$\theta_1$と$eta_2$は直交し、$\theta_2$と$eta_1$は直交する。 Because two coordinate systems are in dual orthogonal relataion, $\theta_1$ and $\eta_2$ are perpendicular and $\theta_2$ and $\eta_1$ are perpendicular.

## 例 Example

- 以下の図は、正規分布の$\theta$座標系格子と$\eta$座標系格子とを、$(m,s)$座標系に描いたものである。The figure below is drawn as its x and y axes are m and s and $\theta$-lattice and $\eta$-lattice are drawn in it.

- 赤線は$\theta_1$一定の点を描いたもの、同様に、緑は$\theta_2$、青は$\eta_1$、水色は$\eta_2$。

```{r,echo=FALSE}
my.ms2theta <- function(m,s){
	theta1 <- m/(s^2)
	theta2 <- -1/(2*s^2)
	return(c(theta1,theta2))
}

my.ms2eta <- function(m,s){
	eta1 <- m
	eta2 <- s^2+m^2
	return(c(eta1,eta2))
}

my.theta2ms <- function(theta1,theta2){
	m <- -theta1/(2*theta2)
	s <- sqrt(-1/(2*theta2))
	return(c(m,s))
}

my.eta2ms <- function(eta1,eta2){
	m <- eta1
	s <- sqrt(eta2-eta1^2)
	return(c(m,s))
}

my.theta2eta <- function(theta1,theta2){
	ms <- my.theta2ms(theta1,theta2)
	my.ms2eta(ms[1],ms[2])
}
my.eta2theta <- function(eta1,eta2){
	ms <- my.eta2ms(eta1,eta2)
	my.ms2theta(ms[1],ms[2])
}
```


```{r}
# 正規分布、(m,s)座標に(theta1,theta2)格子
# 
theta1 <- seq(from=1/2,to=1000,by=1/2)
theta2 <- -seq(from=1/2,to=1000,by=1/2)
eta1 <- seq(from=-1,to=1,by=1/(2^4))
eta2 <- seq(from=0,to=2,by=1/(2^4))

fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)
}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}
```

- 描かれた線はすべて「直線」である

- $\theta_1$一定の線()と$\eta_2$一定の線(水色)は直交し、$\theta_2$一定の線()と$\eta_1$一定の線()も直交している

# Kullback-Leibler ダイバージェンス Kullback-Leibler divergence

## 曲がっているので距離は定まらない Due to non-flatness, no distance is difined

- 双対平坦ではあるが、単純に平坦ではないので、いわゆる距離は定められない。Although dually flat, the space of distributions is not "simply flat". Therefore "distance" in a regular sense can not be defined. 

## ダイバージェンス Divergence

- ダイバージェンスと呼ばれる、距離と似た性質を持つものが定義できる。 Divergence that shares some features with distance can be defined.

- 双対平坦座標系にそのようなものを定めると、Kulback-Leibler divergence (KL divergence) となる。The divergence for dually flat coordinate systems is called Kullback-Leibler divergence (KL divergence).

- 2点 P,Qの確率分布がp,qであるとすると、PからQへのKLdは以下の式で表される。Assume two points P and Q, corresponding to distributions p and q, respectively. KLd from P to Q is expressed as;

$$
KLd(P \to Q) = \int p \frac{\log{p}}{\log{q}} dx
$$

- 式の非対称性からわかるように、QからPに変えると値が異なる。The formula is asymmetric, that means the exchange of P and Q changes the value of divergence.

$$
KLd(Q \to P) = \int p \frac{\log{p}}{\log{q}} dx \ne KLd(P \to Q)
$$

## 例 Example

- 2つの正規分布P(m=0.4,s=0.4)とQ(m=0.6,0.8)を考える。 Assume two normal distributions P(m=0.4,s=0.4) and Q(m=0.6,0.8).

```{r}
x <- seq(from=-10,to=10,length=10^4)
p <- dnorm(x,0.4,0.4)
q <- dnorm(x,0.6,0.8)

matplot(x,cbind(p,q),type="l")
```

- KLd の計算 Calculation of KLd.

```{r}
KLdPQ <- sum(p*(log(p)-log(q))) * (x[2]-x[1])
KLdQP <- sum(q*(log(q)-log(p))) * (x[2]-x[1])
KLdPQ
KLdQP
```
## 正規分布の2次元の広がり Draw normal distributions in 2d space.
$$
\theta_1 = \frac{m}{s^2}\\
\theta_2 = -\frac{1}{2s^2}\\
m = -\frac{\theta_1}{2\theta_2}\\
s^2 = -\frac{1}{2\theta_2}
$$


$$
\eta_1 = m\\
\eta_2 = m^2+s^2\\
m = \eta_1\\
s^2 = \eta_2- \eta_1^2
$$

- 点$P =(m_p=0.4,s_p=0.4)$を取る Point $P =(m_p=0.4,s_p=0.4)$ .

$$
\theta_p1 = \frac{m_p}{s_p^2} = \frac{0.4}{0.16}=2.5\\
\theta_p2 = -\frac{1}{2 \times s_p^2} = -\frac{1}{2\times 0.16}\\
\eta_p1 = m_p = 0.4\\
\eta_p2 = m_p^2 + s_p^2 = 0.32
$$

- 点Pが乗っている赤い線は$\theta_1=2.5$であり、点の近くの緑の線は$\theta_2 = -6$、2本の青い線($m=\eta_1=0.3750,0.4375$)の間に位置し、水色の円では半径$0.3125$の付近にある。 A red line in whic P is, is $\theta_1=2.5$. The green line near P is $\theta_2 = -6$, two blue lines sandwiching P are $m=\eta_1=0.3750$ and $0.4375$. The light blue arc close to P is a part of a circle with radius 0.3125.

```{r}
# 正規分布、(m,s)座標に(theta1,theta2)格子
# 
theta1 <- seq(from=1/2,to=1000,by=1/2)
theta2 <- -seq(from=1/2,to=1000,by=1/2)
eta1 <- seq(from=-1,to=1,by=1/(2^4))
eta2 <- seq(from=0,to=2,by=1/(2^4))

fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
P.theta <- my.ms2theta(P.ms[1],P.ms[2])
points(P.ms[1],P.ms[2],pch=20,cex=3)
text(P.ms[1]+0.1,P.ms[2],"P")
```

- 第2点$Q =(m_p=0.6,s_p=0.8)$を取る The 2nd point $Q =(m_p=0.6,s_p=0.8)$
$$
\theta_q1 = \frac{m_q}{s_q^2} = \frac{0.6}{0.64}=0.9375\\
\theta_q2 = -\frac{1}{2\times s_q^2} = -\frac{1}{2 \times 0.64}\\
\eta_q1 = m_p = 0.6\\
\eta_q2 = m_p^2 + s_p^2 = 1
$$

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=3)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=3,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
```

- 点Pから点Qへと向かうにあたり、$\theta_1$が一定()か$\theta_2$()が一定か、$\eta_1$()が一定か$\eta_2$(水色)が一定な『直線』に沿って進むことにする。Assume you move from P to Q. You take a path on which $\theta_1$ or $\theta_2$ or $\eta_1$ or $\eta_2$ is constant, which is red, green, blue or light blue. This means you move "straight".

- このとき、$\theta_1$一定の線()と$\eta_2$一定の線(水色)は直交しているから、そのどちらかが一定になるように進む方法と Because red lines and light blue lines ($\theta_1$ constant and $\eta_2$ constant) are perpendicular, you can move only using these.

- もう片方の直交ペア$\theta_2$一定の線()と$\eta_1$一定の線()のみを使って進む方法の2つを考える。Another way is to use only greens and blues, that are $\theta_2$ constant or $\eta_1$ constant.

- この2通りの進み方だと、2点間を結ぶ線を斜辺とする「直角三角形」の直角を挟む2辺となるからこのようにする。そして、そのようにするのは、斜辺の長さの取り方として好都合であるからである。どのように好都合かは後述する。The reason why we take these combinations is because we want to have right-angled triangles that are used later. 

- $\theta_1$, $\eta_2$の場合 Case of $\theta_1$ or $\eta_2$

- 太い黒い線に沿って「長方形」の2辺を取ると、片方は$\theta$座標系の直線、もう片方は$\eta$座標系の直線。Thick black lines make "rectangle", the edges of which are straight in $\theta$ or $\eta$ coordinate systems.

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
points(t,sqrt(t/P.thetas[1]),type="l",col=1,lwd = 3)
points(t,sqrt(t/Q.thetas[1]),type="l",col=1,lwd = 3)

points(sqrt(P.etas[2]) * cos(t),sqrt(P.etas[2])* sin(t),type="l",col=1,lwd = 3)
points(sqrt(Q.etas[2]) * cos(t),sqrt(Q.etas[2])* sin(t),type="l",col=1,lwd = 3)
```


- 2本の直線の交点が2つある。$R_1,R_2$とすれば、その座標は Let $R_1$ and $R_2$ denote two crossing points.

$$
R_1(\theta_1=\theta_1^P=2.5,\eta_2=\eta_2^Q = 1)\\
R_2(\theta_1 =\theta_1^Q = \frac{0.6}{0.64},\eta_2=\eta_2^P = 0.32)
$$

- $\theta_1$と$\eta_2$とでできた座標なので、「混合座標」と呼ばれる。この2点の$m,s$を出しておこう。  Their coordinates can be expressed as mixture of $\theta$ and $\eta$; this is why these coordinates are called "mixture coordinates".

```{r}
my.mix2ms12 <- function(theta1,eta2){
  s.sq <- (-1 + sqrt(1+4*theta1^2*eta2))/(2*theta1^2)
  s <- sqrt(s.sq)
  m <- theta1*s.sq
  return(c(m,s))
}
my.mix2ms21 <- function(eta1,theta2){
  m <- eta1
  s.sq <- -1/(2*theta2)
  s <- sqrt(s.sq)
  return(c(m,s))
}
R1.ms <- my.mix2ms12(P.thetas[1],Q.etas[2])
R2.ms <- my.mix2ms12(Q.thetas[1],P.etas[2])
R1.ms
R2.ms
```

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
points(t,sqrt(t/P.thetas[1]),type="l",col=1,lwd = 3)
points(t,sqrt(t/Q.thetas[1]),type="l",col=1,lwd = 3)

points(sqrt(P.etas[2]) * cos(t),sqrt(P.etas[2])* sin(t),type="l",col=1,lwd = 3)
points(sqrt(Q.etas[2]) * cos(t),sqrt(Q.etas[2])* sin(t),type="l",col=1,lwd = 3)

points(R1.ms[1],R1.ms[2],pch=20,cex=5,col=3)
points(R2.ms[1],R2.ms[2],pch=20,cex=5,col=4)
text(R1.ms[1]+0.1,R1.ms[2],"R1")
text(R2.ms[1]-0.1,R2.ms[2],"R2")
```



- $\theta_2$, $\eta_1$の場合 Another case with $\theta_2$ or $\eta_1$.

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
abline(h=sqrt(-1/(2*P.thetas[2])),col=1,lwd=3)
abline(h=sqrt(-1/(2*Q.thetas[2])),col=1,lwd=3)

abline(v=P.etas[1],col=1,lwd=3)
abline(v=Q.etas[1],col=1,lwd=3)

```

- 2本の直線の交点が2つある。$R_3,R_4$とすれば、その座標は Let $R_3,R_4$ denote the crossing points.

$$
R_3(\eta_1=\eta_1^P=0.4,\theta_2=\theta_2^Q = -\frac{1}{2\times0.64})\\
R_4(\eta_1 =\eta_1^Q = 0.6,\theta_2=\theta_2^P = -\frac{1}{2\times 0.16})
$$

- $\theta_2$と$\eta_1$とでできた座標なので、これも「混合座標」と呼ばれる。この2点の$m,s$を出しておこう。 Again they are mixture.

```{r}

R3.ms <- my.mix2ms21(P.etas[1],Q.thetas[2])
R4.ms <- my.mix2ms21(Q.etas[1],P.thetas[2])
R3.ms
R4.ms
```

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
abline(h=sqrt(-1/(2*P.thetas[2])),col=1,lwd=3)
abline(h=sqrt(-1/(2*Q.thetas[2])),col=1,lwd=3)

abline(v=P.etas[1],col=1,lwd=3)
abline(v=Q.etas[1],col=1,lwd=3)

points(R3.ms[1],R3.ms[2],pch=20,cex=5,col=5)
points(R4.ms[1],R4.ms[2],pch=20,cex=5,col=6)
text(R3.ms[1]-0.1,R3.ms[2],"R3")
text(R4.ms[1]+0.1,P.ms[2],"R4")
```

- $\theta_1$,$\eta_2$のペアでの「長方形」も$\theta_2$,$\eta_1$のペアでの「長方形」も、次の意味で「同じ」である。The 1st pair $\theta_1,\eta_2$ generated a rectangle and $\theta_2,\eta_1$ pair also generated a rectangle.




```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
points(t,sqrt(t/P.thetas[1]),type="l",col=4,lwd = 3)
points(t,sqrt(t/Q.thetas[1]),type="l",col=4,lwd = 3)

points(sqrt(P.etas[2]) * cos(t),sqrt(P.etas[2])* sin(t),type="l",col=1,lwd = 3)
points(sqrt(Q.etas[2]) * cos(t),sqrt(Q.etas[2])* sin(t),type="l",col=1,lwd = 3)


abline(h=sqrt(-1/(2*P.thetas[2])),col=4,lwd=3)
abline(h=sqrt(-1/(2*Q.thetas[2])),col=4,lwd=3)

abline(v=P.etas[1],col=1,lwd=3)
abline(v=Q.etas[1],col=1,lwd=3)

points(R1.ms[1],R1.ms[2],pch=20,cex=5,col=3)
points(R2.ms[1],R2.ms[2],pch=20,cex=5,col=4)
points(R3.ms[1],R3.ms[2],pch=20,cex=5,col=5)
points(R4.ms[1],R4.ms[2],pch=20,cex=5,col=6)
text(R1.ms[1]+0.1,R1.ms[2],"R1")
text(R2.ms[1]-0.1,R2.ms[2],"R2")
text(R3.ms[1]-0.1,R3.ms[2],"R3")
text(R4.ms[1]+0.1,R4.ms[2],"R4")
```

## 2点間のダイバージェンス Divergence between two points.

- 2点間の距離のようなものであるダイバージェンスを考える。2点間に直線のようなものを取って、その長さが知りたい。Divergence that something like asquared distance is introduced as a measure between two points. Distance is a measure along straight lines and divergence should be also defined along "straight lines".

- 距離の二乗に相当する量である。Divergence is like square of distance. 

- この空間には、$\theta$座標系での直線も引けるし、$\eta$座標系での直線も引ける。We have two kinds of straight lines in $\theta$ and $\eta$ coordinates systems.

- それぞれの座標系の直線は、線形に表せるから、2点$P=(\theta_1^P,\theta_2^P)$と$Q=(\theta_1^Q,\theta_2^Q)$との間の$\theta$座標系な直線は以下のように表せる。Straignt lines are expressed as linear sum of each coordinates and the coordinates of lines in $\theta$ system are:

$$
(\theta_1,\theta_2) = (\theta_1^P,\theta_2^P) + t (\theta_1^Q-\theta_1^P, \theta_2^Q-\theta_2^P)
$$

- また、$\eta$座標系の直線は次のように表せる。Same for straight lines in $\eta$ system.

$$
(\eta_1,\eta_2) = (\eta_1^P,\eta_2^P) + t (\eta_1^Q-\eta_1^P, \eta_2^Q-\eta_2^P)
$$

- 二通りの「直線」を描く(青太線、赤破線)。Two straight lines are drawn with thick blue line and red broken line, which are ovelapped.

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
points(t,sqrt(t/P.thetas[1]),type="l",col=4,lwd = 3)
points(t,sqrt(t/Q.thetas[1]),type="l",col=4,lwd = 3)

points(sqrt(P.etas[2]) * cos(t),sqrt(P.etas[2])* sin(t),type="l",col=1,lwd = 3)
points(sqrt(Q.etas[2]) * cos(t),sqrt(Q.etas[2])* sin(t),type="l",col=1,lwd = 3)


abline(h=sqrt(-1/(2*P.thetas[2])),col=4,lwd=3)
abline(h=sqrt(-1/(2*Q.thetas[2])),col=4,lwd=3)

abline(v=P.etas[1],col=1,lwd=3)
abline(v=Q.etas[1],col=1,lwd=3)

points(R1.ms[1],R1.ms[2],pch=20,cex=5,col=3)
points(R2.ms[1],R2.ms[2],pch=20,cex=5,col=4)
points(R3.ms[1],R3.ms[2],pch=20,cex=5,col=5)
points(R4.ms[1],R4.ms[2],pch=20,cex=5,col=6)
text(R1.ms[1]+0.1,R1.ms[2],"R1")
text(R2.ms[1]-0.1,R2.ms[2],"R2")
text(R3.ms[1]-0.1,R3.ms[2],"R3")
text(R4.ms[1]+0.1,R4.ms[2],"R4")

t <- seq(from=0,to=1,length=100)
line.theta1 <- P.thetas[1] + t * (Q.thetas[1]-P.thetas[1])
line.theta2 <- P.thetas[2] + t * (Q.thetas[2]-P.thetas[2])
line.eta1 <- P.etas[1] + t * (Q.etas[1]-P.etas[1])
line.eta2 <- P.etas[2] + t * (Q.etas[2]-P.etas[2])
line.theta.ms <- matrix(0,length(t),2)
line.eta.ms <- matrix(0,length(t),2)
for(i in 1:length(t)){
  line.theta.ms[i,] <- my.theta2ms(line.theta1[i],line.theta2[i])
  line.eta.ms[i,] <- my.eta2ms(line.eta1[i],line.eta2[i])
}
points(line.theta.ms[,1],line.theta.ms[,2],type="l",col=4,lwd=5)
points(line.eta.ms[,1],line.eta.ms[,2],type="l",col=2,lwd=2,lty=2)
```

## 二通りのダイバージェンスとKL divergence, Relation between KL divergence and two types of divergences

- 2点P,Qの間のダイバージェンスには、次のような2通りの値がある。PとQとを交換した形になっていて、"$P||Q$"の、もしくは"$Q||P$"のKL divergenceと呼ばれる。Between two points P and Q, two divergence values are defined. KL divergence with "$P||Q$" and "$Q||P$".

$$
KLd(P||Q)=\int p (\log(p)-\log(q))dx\\
KLd(Q||P) = \int q (\log(q)-\log(p))dx
$$

- 情報幾何では、これに次のような4つの意味を持たせる。Information geometry interprets these two in four ways.
  - $D^{\theta}(P\to Q)$:PからQへ、$\theta$的に($\alpha=1$的に)取ったダイバージェンス Divergence from P to Q in terms of $\theta, or $\alpha=1$.
  - $D^{\eta}(P\to Q)$:PからQへ、$\eta$的に($\alpha=-1$的に)取ったダイバージェンス Divergence from P to Q in terms of $\eta, or $\alpha=-1$.
  - $D^{\theta}(Q\to P)$:QからPへ、$\theta$的に($\alpha=1$的に)取ったダイバージェンス Divergence from Q to P in terms of $\theta, or $\alpha=1$.
  - $D^{\eta}(Q\to P)$:QからPへ、$\eta$的に($\alpha=-1$的に)取ったダイバージェンス Divergence from Q to P in terms of $\eta, or $\alpha=-1$.

- これが、以下のような対応関係を持っている。See following relation among them.

$$
D^{\eta}(P\to Q) = D^{\theta}(Q \to P) = KLd(P||Q)\\
D^{\theta}(P\to Q) = D^{\eta}(Q \to P) = KLd(Q||P)\\
$$

- ダイバージェンスを計算してみる。Calculation of divergence values.

```{r}
my.divergence <- function(dP,dQ,dx,alpha=-1){
  if(alpha==-1){
    ret <- sum(dP * (log(dP)-log(dQ))) * dx
  }else{
    ret <- sum(dQ * (log(dQ)-log(dP))) * dx
  }
  return(ret)
}
x <- seq(from=-10,to=10,length=10^4)

dP <- dnorm(x,P.ms[1],P.ms[2])
dQ <- dnorm(x,Q.ms[1],Q.ms[2])

dx <- x[2]-x[1]
# P->Q alpha=1
my.divergence(dP,dQ,dx,alpha=1)
# P->Q alpha=-1
my.divergence(dP,dQ,dx,alpha=-1)
# Q->P alpha=1
my.divergence(dQ,dP,dx,alpha=1)
# Q->P alpha=-1
my.divergence(dQ,dP,dx,alpha=-1)
```

## 拡張ピタゴラスの定理 Extended Pythagorean theorem

- $\theta$,$\eta$は双直交なので、ピタゴラスの定理のような性質が、上記のダイバージェンスには存在する。 Because $\theta$ and $\eta$ are dual orthogonal, divergeces for them have a feature related with Pythagorean theorem.

- ピタゴラスの定理は、直角三角形の3辺の長さに$a^2+b^2=c^2$という関係があることを言う。Pythagorean theorem shows relation among length of three edges of right-angled triangle, $a^2+b^2=c^2$.

- 双対座標系での拡張ピタゴラスの定理は次のような関係のことを言う。The extended Pythagorean theorem in dual coordinate systems is as follow.

- $P$と$R_i$とは$\theta$での「直線(測地線)」で結ばれ、$R_i$と$Q$とは$\eta$での「直線(測地線)」で結ばれているとき、以下の関係にある。具体的には、$R_1,R_4$がこれに相当する。Assume a point $R_i$ for P and Q, where $R_i$ should meet; the "straight"" line between P and $R_i$ is $\theta$' geodesics and the "straight" line between Q and $R_i$ is $\eta$'s geodesics, for example $R_1$ and $R_4$, then;
$$
D^{\theta}(P\to Q) = D^{\theta}(P\to R_i) + D^{\theta}(R_i\to Q)\\
$$


```{r}
dR1 <- dnorm(x,R1.ms[1],R1.ms[2])
dR2 <- dnorm(x,R2.ms[1],R2.ms[2])
dR3 <- dnorm(x,R3.ms[1],R3.ms[2])
dR4 <- dnorm(x,R4.ms[1],R4.ms[2])

my.divergence(dP,dQ,dx,alpha=1)

my.divergence(dP,dR1,dx,alpha=1) + my.divergence(dR1,dQ,dx,alpha=1)
my.divergence(dP,dR4,dx,alpha=1) + my.divergence(dR4,dQ,dx,alpha=1)
```

- $P$と$R_i$とは$\eta$での「直線(測地線)」で結ばれ、$R_i$と$Q$とは$\theta$での「直線(測地線)」で結ばれているとき、以下の関係にある。具体的には、$R2,R3$がこれに相当する。Assume a point $R_i$ for P and Q, where $R_i$ should meet; the "straight"" line between P and $R_i$ is $\eta$' geodesics and the "straight" line between Q and $R_i$ is $\theta$'s geodesics, for example $R_2$ and $R_3$, then;
$$
D^{\eta}(P\to Q) = D^{\eta}(P\to R_i) + D^{\eta}(R_1\to Q)\\
$$


```{r}
my.divergence(dP,dQ,dx,alpha=-1)

my.divergence(dP,dR2,dx,alpha=-1) + my.divergence(dR2,dQ,dx,alpha=-1)
my.divergence(dP,dR3,dx,alpha=-1) + my.divergence(dR3,dQ,dx,alpha=-1)
```

- 斜辺を共有する直角三角形が円周(球面)をなすことを利用すると、双対平坦空間では、2つの円周(球面)が作られる。Right-angled triangles sharing a oblique segment make a circle(sphere surface). Therefore in the dual flat space, two circles(spheres) appear.

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,xlim=c(0,0.6),ylim=c(0.4,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=4)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=5)
}

P.ms <- c(0.4,0.4)
points(P.ms[1],P.ms[2],pch=20,cex=5)
text(P.ms[1]+0.1,P.ms[2],"P")
Q.ms <- c(0.6,0.8)
points(Q.ms[1],Q.ms[2],pch=20,cex=5,col=2)
text(Q.ms[1]+0.1,Q.ms[2],"Q")
P.thetas <- my.ms2theta(P.ms[1],P.ms[2])
P.etas <- my.ms2eta(P.ms[1],P.ms[2])
Q.thetas <- my.ms2theta(Q.ms[1],Q.ms[2])
Q.etas <- my.ms2eta(Q.ms[1],Q.ms[2])
points(t,sqrt(t/P.thetas[1]),type="l",col=4,lwd = 3)
points(t,sqrt(t/Q.thetas[1]),type="l",col=4,lwd = 3)

points(sqrt(P.etas[2]) * cos(t),sqrt(P.etas[2])* sin(t),type="l",col=1,lwd = 3)
points(sqrt(Q.etas[2]) * cos(t),sqrt(Q.etas[2])* sin(t),type="l",col=1,lwd = 3)


abline(h=sqrt(-1/(2*P.thetas[2])),col=4,lwd=3)
abline(h=sqrt(-1/(2*Q.thetas[2])),col=4,lwd=3)

abline(v=P.etas[1],col=1,lwd=3)
abline(v=Q.etas[1],col=1,lwd=3)

points(R1.ms[1],R1.ms[2],pch=20,cex=5,col=3)
points(R2.ms[1],R2.ms[2],pch=20,cex=5,col=4)
points(R3.ms[1],R3.ms[2],pch=20,cex=5,col=5)
points(R4.ms[1],R4.ms[2],pch=20,cex=5,col=6)
text(R1.ms[1]+0.1,R1.ms[2],"R1")
text(R2.ms[1]-0.1,R2.ms[2],"R2")
text(R3.ms[1]-0.1,R3.ms[2],"R3")
text(R4.ms[1]+0.1,R4.ms[2],"R4")

t <- seq(from=0,to=1,length=100)
line.theta1 <- P.thetas[1] + t * (Q.thetas[1]-P.thetas[1])
line.theta2 <- P.thetas[2] + t * (Q.thetas[2]-P.thetas[2])
line.eta1 <- P.etas[1] + t * (Q.etas[1]-P.etas[1])
line.eta2 <- P.etas[2] + t * (Q.etas[2]-P.etas[2])
line.theta.ms <- matrix(0,length(t),2)
line.eta.ms <- matrix(0,length(t),2)
for(i in 1:length(t)){
  line.theta.ms[i,] <- my.theta2ms(line.theta1[i],line.theta2[i])
  line.eta.ms[i,] <- my.eta2ms(line.eta1[i],line.eta2[i])
}
points(line.theta.ms[,1],line.theta.ms[,2],type="l",col=4,lwd=5)
points(line.eta.ms[,1],line.eta.ms[,2],type="l",col=2,lwd=2,lty=2)

KLD1 <- my.divergence(dP,dQ,dx,alpha=1)
KLD2 <- my.divergence(dP,dQ,dx,alpha=-1)

dx <- x[2]-x[1]
n.r <- 10000
r.ms <- cbind(runif(n.r)*1,runif(n.r)*0.7+0.3)

kld1 <- kld2 <- rep(0,length(r.ms[,1]))

for(i in 1:length(r.ms[,1])){
	tmp.d <- dnorm(x,r.ms[i,1],r.ms[i,2])
	kld1[i] <- my.divergence(dP,tmp.d,dx,alpha=1) + my.divergence(tmp.d,dQ,dx,alpha=1)

	kld2[i] <- my.divergence(dP,tmp.d,dx,alpha=-1) + my.divergence(tmp.d,dQ,dx,alpha=-1)
	col <- 1
	#if(col>1)col<-1
	#if(col<0)col<-0
	if(abs(kld1[i]-KLD1)<0.01)col <- 2
	if(abs(kld2[i]-KLD2)<0.01)col <- 3
	if(col>1)points(r.ms[i,1],r.ms[i,2],pch=20,col=col)
}


```

## 練習問題 Exercise

- PとQとを入れ替えて、拡張ピタゴラスの定理を確かめてみよ。Exchange P and Q and check the theorem.

- 上記で選ばれた$R_i$,$\alpha$の組み合わせと異なる組み合わせでは、拡張ピタゴラスの定理が成立しないことを確かめよ。Check the sum does not fit if you apply combinations of $R_i$ and $\alpha$ not mentioned above.

# 推定 Estimation

- 今、観測データが指し示す確率分布が点$P$に相当しているとする。Assume your data set indicates a distribution corresponding to a point P.

- モデルを考えていてそのモデルでは、$\eta_2=a$という固定値だとする。そのときに点Pから、$\eta_2=a$という部分空間への最短ダイバージェンスをとるとする。You want to apply a model where $\eta_2$ value is fixed at a value $a$. You want to fit your data set or $P$ to a subspace where $\eta_2=a$. You can draw a shortest path from $P$ to the subspace. The shortest path should have "shortest" divergence.

- そのような点は、$R_1$である。そのような点R1を探すにあたって、$\theta_1$が一定の直線に沿って探せばよい。In case of the expample shown above, the point to be estimated is $R_1$. You can reach to $R_1$ by taking a straight line of $\theta_1$ being constant.

- これは推定の情報幾何的解釈($\theta$,$\eta$座標系を用いた解釈)の一例である。This is an interpretation of estimation taski in information geometry world with $\theta$ and $\eta$ coordinate systems.

2017-03-08 統計遺伝学分野で生き抜くための数学〜精選18

[][][]指数型分布族と双対座標系

---
title: "指数型分布族と双対座標系 Exponential Family and dual coordinate systems"
author: "ryamada"
date: "2017年3月6日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# はじめに Introduction

- 正規分布を例に都合のよい2つの座標系があることを示し、それらは平坦で、双直交であると説明した。

$$
\theta_1 = \frac{m}{s^2}\\
\theta_2 = -\frac{1}{2s^2}\\
m = -\frac{\theta_1}{2\theta_2}\\
s^2 = -\frac{1}{2\theta_2}
$$


$$
\eta_1 = m\\
\eta_2 = m^2+s^2\\
m = \eta_1\\
s^2 = \eta_2- \eta_1^2
$$

- この座標系のとり方を説明するために、分布の指数型表現というものを使う。

# 正規分布の指数型表現

$$
P(x|m,s) = \frac{1}{\sqrt{2\pi s^2}}e^{-\frac{(x-m)^2}{2s^2}}
$$
- 対数を取り、$x$とパラメタが関係する部分と、パラメタのみにかかわる部分(と$x$のみにかかわる部分(以下の例ではない))とに分ける。



$$
\log{P(x|m,s)} = - \frac{(x-m)^2}{2s^2} -\frac{1}{2}\log{2\pi s^2}\\
= \frac{m}{s^2} x - \frac{1}{2s^2}x^2 - (\frac{m^2}{2s^2} + \frac{1}{2}\log{2\pi s^2})\\
= \begin{pmatrix}\frac{m}{s^2},- \frac{1}{2s^2} \end{pmatrix} \cdot \begin{pmatrix} x \\ x^2 \end{pmatrix} -  (\frac{m^2}{2s^2} + \frac{1}{2}\log{2\pi s^2})
$$

- 次のように書き表すことにする。

$$
\log{P(x|\theta_1,\theta_2)} = \begin{pmatrix}\theta_1,\theta_2 \end{pmatrix} \cdot \begin{pmatrix} x \\ x^2 \end{pmatrix} -  A(\theta_1,\theta_2)\\
= \begin{pmatrix}\theta_1,\theta_2 \end{pmatrix} \cdot \begin{pmatrix} f_1(x) \\ f_2(x) \end{pmatrix} -  A(\theta_1,\theta_2)\\
\theta_1 = \frac{m}{s^2}\\
\theta_2 = -\frac{1}{2s^2}\\
A(\theta_1,\theta_2) = -\frac{\theta_1^2}{4\theta_2} +\frac{1}{2}\log{\frac{\pi}{\theta_2}}\\
f_1(x) = x\\
f_2(x) = x^2
$$

- 以下のように、指数関数の形で表されるので、「指数型」と呼ばれる。

$$
P(x|\theta_1,\theta_2) = e^{\begin{pmatrix}\theta_1,\theta_2 \end{pmatrix} \cdot \begin{pmatrix} f_1(x) \\ f_2(x) \end{pmatrix} -  A(\theta_1,\theta_2)}
$$

- 平坦で双対座標系の片方となる$\theta$座標系は、分布関数を指数型で表したときに、$x$の関数($f_i(x))$の係数であることが知られている

- 他方、$\eta$座標系は、$\theta$座標系と双対な関係にある座標系であるが、どのようなものがそれに相当するかというと、$x$の関数として現れた$f_i(x)$の期待値であることが知られている。

- 正規分布の場合には$f_1(x)=x$と$f_2(x)=x^2$のそれぞれの期待値である。平均$m$、標準偏差$s$の正規分布の$x$の期待値は$m$そのものであるし、$(x-m)^2$の期待値が$s^2$であるから、$x^2$の期待値は、$m^2+s^2$である。

# $\eta$座標系

- 指数型表現をしたときの$x$の関数の係数が$\theta$座標であり、$x$の関数の期待値が$\eta$座標である。

- $\eta$座標は次のようにも表せる。分布関数の$x$によらない成分、$\theta$のみによる成分に「双対対応」するのが$\eta$であり、その「双対対応」というのは偏微分をとることである、というように読める。

$$
\eta_i = \frac{\partial A(\theta)}{\partial \theta_i}
$$
- これは、対数尤度関数を$\theta_i$で偏微分したものの期待値が0であることを利用することで、以下のように示せる。

$$
\frac{\partial \log{p}}{\partial \theta_i} = f_i(x) - \frac{\partial A(\theta)}{\partial \theta_i}\\
\int \frac{\partial \log{p}}{\partial \theta_i} \times p dx = \int (f_i(x) - \frac{\partial A(\theta)}{\partial \theta_i} ) \times pdx\\
\int \frac{\partial \log{p}}{\partial \theta_i} \times p dx=\int f_i(x) \times p dx - \frac{\partial A(\theta)}{\partial \theta_i} \int p dx\\
\int \frac{1}{p}\frac{\partial p}{\partial \theta_i} \times p dx = E[f_i(x)] - \frac{\partial A(\theta)}{\partial \theta_i} \times 1\\
\int \frac{\partial p}{\partial \theta_i} dx = E[f_i(x)] - \frac{\partial A(\theta)}{\partial \theta_i}\\
0 = E[f_i(x)] - \frac{\partial A(\theta)}{\partial \theta_i}
$$
## 練習問題

### 指数型分布族の例 Examples of exponential family

非常に多くの理論確率分布が含まれる。どのような分布が含まれるか確認せよ。https://en.wikipedia.org/wiki/Exponential_family 

There are very many theoretical probability functions in exponental family. See the URL above.

### Canonical formの導出 Derive canonical form

2つの分布をhttps://en.wikipedia.org/wiki/Exponential_family より選び、その指数型表現を通常のパラメタ表現から導け。

Select two distributions from the site and derive their exponential form.

### ヤコビ行列 Jacobian matrix

正規分布のパラメタ変換におけるヤコビ行列を求めよ。 We changed the parameters of normal distribution. Show its Jacobian matrix. 
$$
\theta_1 = \frac{\mu}{\sigma^2}\\
\theta_2 = - \frac{1}{2\sigma^2}
$$

2017-03-07 統計遺伝学分野で生き抜くための数学〜精選17

[][][][]双対平坦、双直交

---
title: "双対平坦・双直交 Dually flat, dually orthogonal"
author: "ryamada"
date: "2017年3月5日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
library(knitr)
opts_chunk$set(echo = TRUE)
library(rgl)
knit_hooks$set(rgl = hook_webgl)
```

# 測地線 Geodesic

- 飛行機で遠くへ、たとえば大阪からロサンゼルスへ、行くとき、両都市の緯度は大体同じなので、その緯線に沿って飛ぶかと言えば、そんなことはなく、北東に向かって飛び始めて、だんだん真東へと向き、南東に向きを変えて到着する。Imagine a flight from Osaka to Los Angels, whose lattitudes are similar. The course is not along the lattitude line but it heads north-east first then heads to east and changes the direction to south-east.


- それが緯線に沿うよりも短距離だから。The flight course is the shortest path.

- 両都市を通る大円に近いコースを取っている。The flight path is along the great circle on which two cities are. 

- この大円コースが球面での「まっすぐな線」である。This great circle course is "the straight line" on the sphere.

- したがって、地球表面に緯度・経度という座標の取り方(パラメタの取り方)をしたときには、「まっすぐな線」というのは、座標系に対しては曲線となることを意味している。From this, we can tell that straight lines are curved in the coordinate system of lattitude and longitude.

- このような曲がった空間でのまっすぐな線が測地線。The straight lines in the curved space are geodesics.

# 平坦 Flatness

- 空間には色々な座標の取り方があるが、ある座標系を取って、その座標をきれいな格子と考えて、地図を描き、その地図において、いわゆる直線で進むことにする。前の例では、緯線・経線を引いて、緯線に沿って進む、というような進み方である。We can assign various coordinate systems to a space. Assume you take paths that are linearly expressed in a coordinate system.

- この進み方が、「たまたま」測地線という意味でも「まっすぐな線」であるとき、この座標系の取り方は、「平坦」である、と言う。If the paths appearing straight in the coordinate system are geodesics, then the coordinate system is flat.

## 平坦とフィッシャー情報量と接続係数 Flatness, Fisher information and connection coefficient

- フィッシャー情報量は、局所の長さ・内積の測り方を決めている。Fisher information defines local elongation and inner product.

- フィッシャー情報量は座標系の取り方に応じて表現される。Fisher information's expression depends on the coordinate systems.

- フィッシャー情報量は局所の長さを決めているので、曲線を引くときにその長さをどのくらいと見積もるかを決める。Fisher information tells local length, therefore it determins the length of curves/lines in the space because the length is integral of local lengths.

- したがって、フィッシャー情報量が(フィッシャー情報量を構成する関数の変化具合が)曲線に沿ってどう変化するかを検討することで、どのような座標系の取り方は、平坦なのかが解る。This means that the changes of Fisher information along a curve have information whether the curves are geodesics and also have information on flatness of coordinate system.

- 実際、平坦ならば、0になるべき要素を接続係数$F_{ki,j}^{(\alpha)}$と呼ぶが、それは以下で示されることからもその様子が見て取れる。Actually connection coeffcicient $F_{ki,j}^{(\alpha)}$ is defined so that $F_{ki,j}^{(\alpha)}$ should be zero when flat. The formula below implicates that connection coefficients are closely related to Fisher information.

$$
g(ij) = \int \frac{\partial 2\sqrt{p}}{\partial \theta_i}\frac{\partial 2\sqrt{p}}{\partial \theta_j} dx\\
= \int  \frac{\partial l^{(\alpha)}}{\partial \theta_i}\frac{\partial l^{(-\alpha)}}{\partial \theta_j} dx\\
$$

$$
F_{ki,j}^{(\alpha)} =\int  \frac{\partial}{\partial \theta_k}\frac{\partial l^{(\alpha)}}{\partial \theta_i}\frac{\partial l^{(-\alpha)}}{\partial \theta_j} dx
$$

- ただし、$l^{(\alpha)}$はフィッシャー情報量をうまく2つの関数の偏微分の積を使って表すために登場した関数。where $l^{(\alpha)}$ stands for function that expresses Fisher information as product of partial derivatives of them.

$$
l^{(\alpha)} = \frac{2}{1-\alpha}p^{\frac{1-\alpha}{2}}; \alpha \ne 1\\
= \log{p}; \alpha = 1
$$


- 接続係数が0になるのは、$\alpha=\pm=1$のときであることが示せる。そして$\alpha=\pm 1$の2つはそれぞれConnection coefficients are 0 only when $\alpha=\pm 1$. 


- そして、確率分布の多くが、このような座標系のペアを持てることが知られているので、そのような確率分布を配した空間・多様体は二重平坦である、と言われる。Almost all probability distributions are known to have a pair of coordinate systems that are both flat in the space of distributions. And the manifolds are called dually flat.

# 二重平坦 Dually flat

## 正規分布の例 Example; Normal distribution

- 平均m、標準偏差sで表現すれば、Formula of normal distribution with mean $m$ and standard deviation $s$.
$$
p(x|m,s) = \frac{1}{\sqrt{2\pi}s}e^{-\frac{(x-m)^2}{2s^2}}
$$
- これに$(m,s)$という二次元座標を与えた、と見ることができる。This formula can be read that 2d coordinates $(m,s)$ are given to the distribution.

- 今、異なる2つの座標系$(\theta_1,\theta_2)$,$(\eta_1,\eta_2)$を与える。We take two different coordinate systems $(\theta_1,\theta_2)$,$(\eta_1,\eta_2)$ as below.

- それぞれ、$\alpha = 1$,$\alpha=-1$に対応する2つの座標系として知られているものである。They are known as coordinate systems corresponding to $\alpha = 1$,$\alpha=-1$.

$$
\theta_1 = \frac{m}{s^2}\\
\theta_2 = -\frac{1}{2s^2}\\
m = -\frac{\theta_1}{2\theta_2}\\
s^2 = -\frac{1}{2\theta_2}
$$

$$
\eta_1 = m\\
\eta_2 = m^2+s^2\\
m = \eta_1\\
s^2 = \eta_2- \eta_1^2
$$

- $(\theta_1,\theta_2)$の格子を描く Let's draw lattice of $(\theta_1,\theta_2)$

- 位置によって、格子が作る伸び縮みの具合は変化しているし、直交しているわけでもない The lattice is deformed and the features of elongation varies among locations and the angles are not perpendicular.

- 今、描図に用いている$(m,s)$が絶対的な座標系ではないので、伸び縮みの変化や角度が「本当のところ」はどうなっているのかは、図からは分からない。それはフィッシャー情報量を各所で調べることによってわかる The coordinate system $(m,s)$ is one of many coodinate systems and it is not the right one. Therefore this drawing does not tell "real" elongation features of $(\theta_1,\theta_2)$ system. If you want to know the "real" elongation, you have to check its Fisher information. 

```{r}
# 正規分布、(m,s)座標に(theta1,theta2)格子
# (theta1,theta2) lattice on the map with (m,s) coordiante system
# s^2 = m/theta1
# s^2 = -1/(2theta2)
theta1 <- seq(from=1/2,to=1000,by=1/2)
theta2 <- -seq(from=1/2,to=1000,by=1/2)
fr <- matrix(c(-1,0,1,1),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,ylim=c(0,2))

t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l")
  points(-t,sqrt(t/theta1[i]),type="l")

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])))
}
abline(h=0)
```

- $(\eta_1,\eta_2)$の格子を描くLet's draw lattice of $(\eta_1,\eta_2)$

- 位置によって、格子が作る伸び縮みの具合は変化しているし、直交しているわけでもない Elongation feaures depends on location. Not perpendicular.

- 今、描図に用いている$(m,s)$が絶対的な座標系ではないので、伸び縮みの変化や角度が「本当のところ」はどうなっているのかは、図からは分からない。それはフィッシャー情報量を各所で調べることによってわかる For the true elongation, check Fisher information.

```{r}
# 正規分布、(m,s)座標に(eta1,eta2)格子
# (eta1,eta2) lattice on the map with (m,s) coordiante system

# m = eta1
# m^2+s^2=eta2
eta1 <- seq(from=-1,to=1,by=1/(2^4))
eta2 <- seq(from=0,to=2,by=1/(2^4))

fr <- matrix(c(-1,0,1,1),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,ylim=c(0,2))

for(i in 1:length(eta1)){
	abline(v=eta1[i])
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l")
}
```

- $(\theta_1,\theta_2)$と$(\eta_1,\eta_2)$の2種類の格子を重ねてみる Draw two coordinate system lattice together.


```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,ylim=c(0,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)

}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=2)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=3)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=3)
}

```

- $\theta_1$と$\eta_1$とを同じ色で、$\theta_2$と$\eta_2$とを同じ色で描きなおしてみる Draw the same but colors are common for $\theta_i$ and $\eta_i$ not for $\theta_1$と$\theta_2$.

```{r}
fr <- matrix(c(-0.5,0,0.5,0.5),byrow=TRUE,2,2)
plot(fr,asp=1,col=0,ylim=c(0,1))
t <- seq(from=0,to=1,length=100) * 2*pi

for(i in 1:length(theta1)){
  points(t,sqrt(t/theta1[i]),type="l",col=2)
  points(-t,sqrt(t/theta1[i]),type="l",col=2)
	#segments(0,0,2,2/theta1[i],col=2)
	#segments(0,0,-2,2/theta1[i],col=2)
}

for(i in 1:length(theta2)){
	abline(h=sqrt(-1/(2*theta2[i])),col=3)
}
for(i in 1:length(eta1)){
	abline(v=eta1[i],col=3)
}
t <- seq(from=0,to=1,length=100) * 2*pi
for(i in 1:length(eta2)){
	points(eta2[i] * cos(t),eta2[i] * sin(t),type="l",col=2)
}

```


# 双対座標系 Dual coordinate system

- 正規分布は2次元多様体であるので、2次元座標を与えられる Coordinate systems of normal distribution is 2d.

- 二つの座標系$(\theta_1,\theta_2)$,$(\eta_1,\eta_2)$を示した Two systems $(\theta_1,\theta_2)$ and $(\eta_1,\eta_2)$ were introduced above.

- この座標系は、正規分布が持つ特別な座標系ペアである They are special for normal distribution.

  - どちらも平坦 Both are flat.
  - 相互に直交(双直交) Theya re mutually orthogonal (dually orthogonal).

- どちらも、平坦であるというのは、座標系を「信じて」まっすぐに進むとそれが測地線になること。このことは、それぞれの座標系が持つ性質であって、2つの座標系の関係を考える必要がない Both are flat, which means that straight paths in either coordinate system are geodesics. This feature is intrinsic for both system and no interaction between two systems.

- 相互に直交(双直交)であることを考えるには、それぞれの座標系の個々の要素に対応付けをする必要がある。$\theta_1$と$\eta_1$とは、ある理由で対応があり、その同じ理由で$\theta_2$と$\eta_2$とも対応がある、と言うような関係である。一般に、k次元の場合も$\theta_i$と$\eta_i$,$i=1,2,...,k$に対応がある The dual orthogonality is the feature between two systems. Therefore we need a piece of rule between the two systems. The elements of $\theta$ system, and the elements of $\eta$ system are in the one-to-one correspondence relation. $\theta_1$ is something specific for $\eta_1$ but not for $\eta_2$.

- このとき、$\theta_i$と$\eta_j$とに関する、偏微分ベクトル(局所の座標方向のベクトルであって、伸び縮みを表すベクトル)同士が、正規直行基底であるかのような性質を持ち、そのことを双直交と言う Now we have mutual relation among $\theta_i$ and $\eta_j$ depending on $i=j$ or $i \ne j$ and the set of partial derivative vectors of $\theta$ system and the set of partial derivative vectors of $\eta$ system are in the relation where the inner products of both system vectors orthonormal. And this relation is called dual orthonormal.

- 正規直行基底であるかのよう、とは、$i =j$のときには、内積が1であり、$i\ne j$のときには内積が0となることである It means that the inner product of the partial derivative vectors of $\theta_i$ and $\eta_i$ is 1 and the inner product of vectors of $\theta_i$ and $\eta_j$, where $i \ne j$, is 0.

- $\theta_i$と$\eta_j$とに関する、偏微分ベクトルと書いたが、
  - $\theta_i$ のそれは$l^{(\alpha=1)}=\log{p}$から来る座標系なので、その偏微分ベクトルは$\frac{\partial \log{p}}{\partial \theta_i}$である
  - 他方$\eta_i$のそれは$l^{(\alpha=-1)}=p$から来る座標系なので、その偏微分ベクトルは$\frac{\partial p}{\partial \eta_i}$。The partial derivative vectors are somewhat different between $\theta$ system and $\eta$ system. For $\theta$ system, differentiate $\log{p}$ and for $\eta$ system, differentiate $p$ itself.
  - 別の言い方もできる。$\int \frac{\partial \log{p}}{\partial \theta_i} \frac{\partial p}{\partial \eta_i} dx$は$\int \frac{\partial \log{p}}{\partial \theta_i} \frac{\partial \log{p}}{\partial \eta_i} \times p dx$と変形できるから、$\log{p}$の$\theta_i$と$\eta_j$とのそれぞれの偏微分の期待値が1または0になる関係でもある Another explanation is also possible. $\int \frac{\partial \log{p}}{\partial \theta_i} \frac{\partial p}{\partial \eta_i} dx$ can be transformed to $\int \frac{\partial \log{p}}{\partial \theta_i} \frac{\partial \log{p}}{\partial \eta_i} \times p dx$. This means $\log{p}$ is partially differentiated with $\theta_i$ and $\eta_j$ and their "expected value" is 1 or 0.
  

  
## 正規分布の例 Example, normal distribution

- 具体例で確認する Confirm what the above descriptions mean with an example.

- ある正規分布$N(m,s^2)$を考え、それを、少し変化させる。変化させる方向は$\theta_1$,$\theta_2$,$\eta_1$,$\eta_2$とすれば、それぞれの微小変化によって正規分布が変わる、$x \in R$全体にわたって増えたり減ったりする Take an instance of normal distribution $N(m,s^2)$ and change its parameters a bit. The direction of change is in $\theta_1$,$\theta_2$,$\eta_1$ or $\eta_2$. With the change in any one of four directions, values of normal distributio for $x \in R$ changes everywhere. 

- その変化の総和(xに関する積分)をパラメタの増分で割ったものが、$\frac{\partial \log{p}}{\partial \theta_i}$だったり、$\frac{\partial p}{\partial \eta_j}$となる The integral of the changes throughout $x$ should divieded by the change in parameter value, that is approxiamtion of derivative.

- それを確かめるために、$(m,s)$と$(\theta_1,\theta_2)$,$(\eta_1,\eta_2)$とを変換する関数や、$(\theta_1,\theta_2)$のときの$\log{N(m,s)}$の値や$(\eta_1,\eta_2)$のときの$N(m,s)$の値を算出する関数を作る First make some utility functions to perform this numeric experiments.

- 座標の相互変換の関数 Functions to convert coordinates among three systems.

```{r}
my.ms2theta <- function(m,s){
	theta1 <- m/(s^2)
	theta2 <- -1/(2*s^2)
	return(c(theta1,theta2))
}

my.ms2eta <- function(m,s){
	eta1 <- m
	eta2 <- s^2+m^2
	return(c(eta1,eta2))
}

my.theta2ms <- function(theta1,theta2){
	m <- -theta1/(2*theta2)
	s <- sqrt(-1/(2*theta2))
	return(c(m,s))
}

my.eta2ms <- function(eta1,eta2){
	m <- eta1
	s <- sqrt(eta2-eta1^2)
	return(c(m,s))
}

my.theta2eta <- function(theta1,theta2){
	ms <- my.theta2ms(theta1,theta2)
	my.ms2eta(ms[1],ms[2])
}
my.eta2theta <- function(eta1,eta2){
	ms <- my.eta2ms(eta1,eta2)
	my.ms2theta(ms[1],ms[2])
}
```

- それぞれの座標系で、$\log{p}$,$p$を返す関数を作る Functions to return $\log{p}$ or $p$ for each coordinate systems.

```{r}
# theta系はlog(p)
my.dnorm.theta <- function(x,theta1,theta2,log=TRUE){
	ms <- my.theta2ms(theta1,theta2)
	dnorm(x,ms[1],ms[2],log=log)
}
# eta系はp
my.dnorm.eta <- function(x,eta1,eta2,log=FALSE){
	ms <- my.eta2ms(eta1,eta2)
	dnorm(x,ms[1],ms[2],log=log)
}
```

- ある特定の正規分布$m0 = 2, s0 = 1$について偏微分ベクトルを求め、双直交性を確認する Take an example normal distribution with $m0 = 2, s0 = 1$ and calculate numeric approximates of partial derivatives for 4 directional chnges.

```{r}
m0 <- 2
s0 <- 1
dx <- 1/1000
x <- seq(from=-10,to=10,by=dx)
p <- dnorm(x,m0,s0)
plot(x,p)
```

- 対応する$\theta,\eta$を求める $\theta,\eta$?
```{r}
thetas <- my.ms2theta(m0,s0)
etas <- my.ms2eta(m0,s0)
```

- $\theta_1$の偏微分を求める。その増減の様子を描く Partial derivative in the direction of  $\theta_1$.
```{r}
d.theta1 <- 0.0001
dp.dtheta1 <- (my.dnorm.theta(x,thetas[1]+d.theta1,thetas[2]) - my.dnorm.theta(x,thetas[1],thetas[2]))/d.theta1

plot(x,dp.dtheta1)
```

- $\theta_2,\eta_1,\eta_2$についても同様に求める Same for $\theta_2,\eta_1,\eta_2$.
```{r}
d.theta2 <- 0.0001

d.eta1 <- 0.0001
d.eta2 <- 0.0001


dp.dtheta2 <- (my.dnorm.theta(x,thetas[1],thetas[2]+d.theta2) - my.dnorm.theta(x,thetas[1],thetas[2]))/d.theta2
plot(x,dp.dtheta2)
dp.deta1 <- (my.dnorm.eta(x,etas[1]+d.eta1,etas[2]) - my.dnorm.eta(x,etas[1],etas[2]))/d.eta1
plot(x,dp.deta1)
dp.deta2 <- (my.dnorm.eta(x,etas[1],etas[2]+d.eta2) - my.dnorm.eta(x,etas[1],etas[2]))/d.eta2
plot(x,dp.deta2)
```

- 双直交性の確認 Dual orthogonality is confirmed as;
```{r}
sum(dp.dtheta1 * dp.deta1) * diff(x)[1]
sum(dp.dtheta1 * dp.deta2) * diff(x)[1]
sum(dp.dtheta2 * dp.deta1) * diff(x)[1]
sum(dp.dtheta2 * dp.deta2) * diff(x)[1]
```