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

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

2016-12-05 指数型分布族の間のKLd

[][][]指数型分布族の間のKLd

  • 一般的な関係式がある
  • 指数分布族の自然パラメタ表現でのポテンシャル関数¥phi(¥theta)によって定まる
  • p(x|¥theta) = e^{C(x)+¥eta(¥theta) T(x) - ¥phi(¥theta)}
  • KL(q||p) = ¥phi(¥theta_p)-¥phi(¥theta_q) - (¥theta_p-¥theta_q)¥cdot ¥nabla ¥phi(¥theta_q)
  • ミスプリ???要、確認

2016-12-04 指数型分布族に近似する

ryamada222016-12-04

[][]指数型分布族に近似する

  • 指数型分布族は便利
  • 情報幾何で便利
  • ただし、複数の正規分布の混合とかが合わないので、不便
  • p(w) = e^{C(w) + F ¥theta - ¥phi(¥theta)}と書ければよいのだから
  • ¥log{p(w)} = C(w) + F ¥theta - ¥phi(¥theta)となるような式を見つければよいのでは…
  • ¥thetaとしてw,w^2,w^3...を取ることにすれば、結局、確率密度関数の対数をとって、それを多項式近似すればよいのでは
  • 無限遠に行けばいくほど確率が0に近づかないといけないので、その点で少し工夫をする

f:id:ryamada22:20161204113321j:image

x <- runif(1000,-10,10)
x <- c(x,c(-20,20))
d <- dnorm(x,2,1) * 0.7 + dnorm(x,4,3) * 0.3

log.d <- log(d)
log.d. <- log.d + rnorm(length(x))*1

log.d.[length(log.d.)] <- log.d.[length(log.d.)-1] <- min(log.d.) - 100


model <- lm(log.d. ~ poly(x,100,raw=TRUE))

op <- par(mar = c(0, 4, 0, 3), oma = c(5, 0, 4, 0), mfcol = c(5,1))

plot(x,log.d.)
plot(x,predict(model))
plot(x,exp(predict(model)))
plot(x,d)
plot(d,exp(predict(model)))

2016-12-03 時系列解析

[]時系列解析

  • メモ
blog <- read.table("blog.txt",col.names=c("統計遺伝学","応用数学"))

ryamada.two <- ts(blog, start=c(2005,1), end=c(2016, 11), frequency=12)

plot(ryamada.two,plot.type="single",col=1:2)


ryamada22 <- ts(blog$統計遺伝学, start=c(2005,1), end=c(2016, 11), frequency=12)
ryamada <- ts(blog$応用数学, start=c(2005,1), end=c(2016, 11), frequency=12)

ryamadaall <- ts(blog$統計遺伝学+blog$応用数学, start=c(2005,1), end=c(2016, 12), frequency=12)

# Seasonal decomposition
fit <- stl(ryamada, s.window="period")
fit22 <- stl(ryamada22, s.window="period")
fitall <- stl(ryamadaall,s.window="period")

ts.data<- ryamada.two
ts.seasonal <- ts(cbind(fit22[[1]][,1],fit[[1]][,1]), start=c(2005,1), end=c(2016, 11), frequency=12)

ts.trend <- ts(cbind(fit22[[1]][,2],fit[[1]][,2]), start=c(2005,1), end=c(2016, 11), frequency=12)

ts.residual <- ts(cbind(fit22[[1]][,3],fit[[1]][,3]), start=c(2005,1), end=c(2016, 11), frequency=12)
op <- par(mfcol = c(2,2))
op <- par(mar = c(0, 4, 0, 3), oma = c(5, 0, 4, 0), mfcol = c(4,1))
plot(ts.data,plot.type="single",col=1:2)
plot(ts.seasonal,plot.type="single",col=1:2)
plot(ts.trend,plot.type="single",col=1:2)
plot(ts.residual,plot.type="single",col=1:2)

plot(fit,main="応用数学",set.pars = NULL)
plot(fit22,main="遺伝統計学",set.pars = NULL)
#plot(fitall)

matplot(cbind("応用数学"=fit[[1]][,1],"統計遺伝学"=fit22[[1]][,1]),type="l",)




# additional plots
monthplot(fit)
monthplot(fit22)
monthplot(fitall)

library(forecast)
seasonplot(fit)
seasonplot(fit22)
seasonplot(fitall)


plot(blog$V1,blog$V2,type="l")

ryamada22.year <- ryamada.year <- rep(0,12)
for(i in 1:12){
	tmp <- (1:12) + (i-1) * 12
	points(blog$V1[tmp],blog$V2[tmp],col=i,pch=20)
	ryamada22.year[i] <- sum(blog$V1[tmp])
	ryamada.year[i] <- sum(blog$V2[tmp])
}

plot(ryamada22.year,ryamada.year)


0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
0	0
1761	50
2401	374
4070	339
3420	345
3314	340
3438	308
5042	501
5004	355
5114	226
4662	0
5283	0
5875	0
6439	0
5908	0
6275	0
5015	0
5198	0
5185	0
5901	0
6045	0
6412	0
5060	0
5391	0
7388	0
9132	0
7705	0
9763	0
7296	0
6216	0
8217	0
9015	0
9604	0
10856	0
8005	0
9175	0
11059	0
10367	0
9370	0
9982	0
7466	0
6063	0
6870	0
7899	0
9560	0
9687	0
6561	0
6446	0
8314	0
7633	0
7185	0
7479	0
6986	0
5830	0
5144	0
6066	0
7262	0
6792	325
6097	1316
6105	1532
8522	2220
8761	3164
7282	2681
8289	3037
7310	2158
5847	1671
7167	1983
8290	2433
8524	2318
8787	2848
5989	2107
6618	2320
7160	2704
7390	2619
7112	2660
8509	3765
8023	3955
6112	3305
7235	3934
9040	4607
9649	4899
9411	5259
7081	4317
6240	4151
8142	5290
8362	5597
8265	5325
8868	5697
7081	4825
5975	4229
7051	5033
8364	5241
7610	5211
7910	6200
6354	4987
6454	4311
8531	5482
9522	6159
9011	5973
9542	6004
7537	4683
6996	4115
7947	5342
8684	5765
9141	5876
9083	6301
7450	5280
7228	5351
8937	5741
8903	6008
8801	6118
9855	7017
7858	5800
6683	5480
9688	5817
8047	5366
8715	5896
8633	5970
6310	4790
7115	4451
7641	5696
7584	5942
7775	6165
8636	5955
7176	6729
7578	5613
7862	5941
9263	6302
9457	6374
10824	7664
8887	6036
8469	5445
7060	5679
9027	6805

2016-11-30 わたしのための情報幾何〜双対平坦

[][][][][]わたしのための情報幾何〜双対平坦

  • 暫定Rmdファイル(ただし、更新されない可能性が大)
---
title: "私のための情報幾何 InformationGeometry4Me"
author: "ryamada"
date: "2016年12月2日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE)
```

# 結論 Conclusion

尤度関数 $p(w,\xi)$ (ただし、$w$は事象、$\xi$はパラメタ(取り方に寄らず、一般的な意味での))には、2つのよいパラメタの取り方$\theta$ と $\eta$がある。

$$
\log{p(w)} = C(w) + F \theta - \phi(\theta)
$$
$$
p(w) = F^{-1}(w) \eta
$$

この二つのパラメタの取り方はには2つの特徴がある。

一つ目の特徴は$\theta$, $\eta$のそれぞれに備わったものであり、もう一つの特徴は、$\theta$と$\eta$との相互関係に関するものである。2つのパラメタの取り方が相互にしていることは、上の2つの式が同一の行列$F$を持っていることから解る。

+ $\theta$, $\eta$のパラメタの取り方をすると、どちらも「平坦」という性質を持つ。「平坦」である、とは、直線が、パラメタの一次線形式で表せるという性質のことである。

+ $\theta$, $\eta$ 座標系は相互に「直交」するという性質を持つ。「直交」するというのは角度が90度である、というよりは、距離を測るときにピタゴラスの定理のようなものを使える、という意味である。その具体的な意味は後述する。

# どのように役に立てるか How it works for me

尤度関数が空間に置かれている(多様体をなしている)。

その多様体上で、個々の尤度関数同士がどれくらい離れているかを測る(尤度比を知るために、最尤推定をするために)りたい。

そのときに、まっすぐな道を通ることが最短距離である性質(平坦)と、ピタゴラスの定理を使えるという性質が便利である。

実際には、$\eta$座標の一部を固定し、残りを自由にすると、すべての$eta$座標が自由な空間から、一部が制約された部分空間への垂線を取るという作業が発生する。

このとき、$\eta$,$\theta$座標系が相互に直交しているので、$\theta$座標の軸に沿った直線が垂線になるという性質があり、垂線とその足を見出すことが容易になる。

# 一目で見る基本概念 Basic concepts @ a glance

尤度関数と情報幾何では、平坦・直線・距離・ピタゴラスの定理を扱うが、ユークリッド幾何的なそれを拡張した上で使っているので、拡張された意味を理解することが必要である。

## 平坦

直線が座標系によって一次関数で表されることを言う。

### 接続
平坦か平坦でないかを定量化したものに「接続」と呼ばれる係数がある。
接ベクトルがある接ベクトル方向に動いたときにまったく変わらなければ、平坦だが、少し変わるとすると、その変わり具合を表現する必要がある。

その変わり具合は色々な方向に起きた変わり具合を併せたものとなる。

結局、$n$個のパラメタで表されている場合には、$n$個の接ベクトルが、$n$個の接ベクトル方向についてどのように変化するかが必要であるから$n^2$ペアの評価が必要となる。

そして、$n^2$のそれぞれは$n$方向成分に分解されるので、$n^2 \times n = n^3$の値の組が、接ベクトルの変化府愛を説明する。

平坦というのは、この$n^3$の係数がすべて0であるような状態のこと、そのようなパラメタの取り方のことである。

### 直線と測地線

平坦な座標系でパラメタの一次線形式で表されるものが直線である。

座標系の取り方を変えると、直線が曲がって見える。

この曲がって見える線が、「最短距離」を表しているときには、「直線」的な意味を持つので、そのような曲線を「測地線」と呼ぶ。


## 距離

### 内積と距離

ユークリッド空間における、2点x=(x_i) と y=(y_i)との距離は

$$
d_{xy}=\sqrt{\sum (x_i-y_i)^2}
$$

と計算される。

$$
d_{xy}^2 = \sum (x_i-y_i)^2 = \sum x_i^2 + \sum y_i^2 - 2 \sum x_i y_i
$$

と書き直し、さらに$\sum x_i y_i$を$x$と$y$との普通の意味での内積であると読めば

$$
d_{xy}^2 = ||x||^2 + ||y||2 - 2 (x,y)
$$

と書ける。さらに$||x||^2$も$x$と$x$自身の内積だと考えれば

$$
d_{xy}^2 = (x,x) + (y,y) - 2(x,y)
$$
となる。

そしてこれを「距離」の定義だと考えることもできる。

普通の意味での内積というのは、ベクトルのペアに関して

$$
(x,y) = x^T G y
$$
と言うベクトルと行列の演算をしているとみなした場合、$G$が単位行列になっていることに相当する。

ユークリッド空間では内積が単位行列で定義されている、と言い換えて、ユークリッド空間でない場合には、内積が単位行列とは異なる正方行列$G$で定義されている、と見ることができる。

この$G$を多様体のリーマン計量と呼ぶ。
また、情報幾何ではフィッシャー情報行列がこのリーマン計量と同じである。

### 距離の拡張としてのKL divergence

距離は内積を使って
$$
d_{xy}^2 = (x,x) + (y,y) - 2(x,y)
$$
と表された。

実際には

ここでの$x$と$y$とは同じ座標系で表されていたのだが、情報幾何では、相互によい関係にある2つの座標系を用いる。

$x$を一つ目の座標系で表し、$y$をもう一つの座標系で表したものとしたときに、一つ目の座標系での$x$と$x$との内積があり、もう一つの座標系で$y$と$y$との内積があり、異なる座標系で表された$x$,$y$とについても$(x,y)$が定義されていれば、上式は計算できる。

実際、

$$
\frac{1}{2} d_{xy}^2 = \frac{1}{2} (x,x) + \frac{1}{2} (y,y) - (x,y) = KLd(x||y)
$$
がKL divergenceと呼ばれる量である。

$(x,x)$,$(y,y)$ は平坦な座標系でのそれなので、「まっすぐ」に測ることができる。

ユークリッド空間で$x$と$y$とが直交していれば、内積$(x,y)$は0になって、

$$
\frac{1}{2} d_{xy}^2 = \frac{1}{2} (x,x) + \frac{1}{2} (y,y)
$$
となるが、これは「ピタゴラスの定理」の性質である。

$\theta$,$\eta$座標系において、両者が直交しているとき、$\theta$座標系での直線$x$と$\eta$座標系での直線$y$とは相互に直交しているから、この式が成り立つ。

それを
$$
KLd(p||r) = KLd(p||q) + KLd(q||r)
$$
と書く。ただし、p-qは$\theta$座標系での直線、q-rは$\eta$座標系での直線である。

## 曲がっていることとリーマン計量と接続

ユークリッド空間では、リーマン計量が単位行列であって、接続係数がすべて0であった。

曲がっている多様体ではリーマン計量が単位行列ではないし、接続係数も0ではない。

逆に言うと、座標系を入れて多様体を考えたり、多様体上の「距離」を測ったりするときには、リーマン計量と接続との両方を考える必要がある。

パラメタの数が$n$個のとき、リーマン計量は$n \times n$行列の形をしているし、接続係数は$n^3$個の値の組になっている。

## 双対

二つの座標系$\theta$,$\eta$が相互にうまい関係になっている話をしている。

この2つは双対とよばれる関係である。お互い様、である。

そのお互い様がどうして、$\log{p(w)}$と$p(w)$と関係するのか、と言ったことを確認するのが、以降の「解説」の記述である。

$$
\log{p(w)} = C(w) + F \theta - \phi(\theta)
$$
$$
p(w) = F^{-1}(w) \eta
$$

# 解説 Details

## 確率密度関数と尤度関数

確率密度関数(確率質量関数)は、パラメタを確率変数の値との関数で

$$
p(w,\xi)
$$
と書ける。

これは尤度関数でもある。

以降は「尤度関数」とよぶことにする。

## 対数尤度関数

$$
L(w,\xi) = \log{p(w,\xi)}
$$
をよく使う。
理由もあるが、それは後述する。

## 尤度関数・対数尤度関数の微分

以下の変換を多用するので、列挙しておく。

$$
L(w,\xi) = \log{p(w,\xi)}
$$
$$
\frac{\partial}{\partial \xi_i} L(w,\xi) = \frac{1}{p(w,\xi)} \frac{\partial}{\partial \xi_i} p(w,\xi)
$$

$$
p(w,\xi) = \frac{\frac{\partial p(w,\xi)}{\partial \xi_i}}{\frac{\partial L(w,\xi)}{\partial \xi_i}}
$$
$$
\frac{\partial}{\partial \xi_i} p(w,\xi) = \frac{\partial L(w,\xi)}{\partial \xi_i} p(w,\xi)
$$


## 対数尤度関数とスコア関数

パラメタ$\xi$の最尤推定は、$\frac{\partial p(w,\xi)}{\partial \xi_i}=0$ を満足する$\xi_i$であるが、$\frac{\partial L(w,\xi)}{\partial \xi_i}=0$ でもよい。

対数尤度関数の方が有益なのは、

$$
E[\frac{\partial L(w,\xi)}{\partial \xi_i}] = 0
$$

であるから。
期待値が0というのは、多数回行ったときに、標本平均を取ることが適切であることを意味するから。

ちなみに$E[\frac{\partial p(w,\xi)}{\partial \xi_i} \ne 0]$ 。

念のため$E[\frac{\partial L(w,\xi)}{\partial \xi_i}] = 0$を確かめる。

$$
E[\frac{\partial L(w,\xi)}{\partial \xi_i}] = \int_\Omega \frac{\partial L(w,\xi)}{\partial \xi_i} p(w,\xi) d\mu = \int_\Omega \frac{1}{p(w,\xi)} \frac{\partial p(w,\xi)}{\partial \xi_i} p(w,\xi) d\mu = \int_\Omega \frac{\partial p(w,\xi)}{\partial \xi_i} d\mu = 0
$$
### スコア関数

このように対数尤度関数の方が尤度関数そのものよりよい性質を持っているので、対数尤度関数の微分には名前がついている。

それがスコア関数である。

$$
\frac{\partial L(w,\xi)}{\partial \xi_i}
$$

## フィッシャー情報量

スコア関数が0のところが最尤推定値である。

最尤推定値が効率よく見つかる条件というのは、スコア関数の値が0の値の周辺で大きく変化していることである。

期待値が0の周りで変化が大きいことは分散が大きいことに相当するので。スコア関数の分散が、0回りで大きいか小さいかを値にすることは有用であると解る。

偏微分をしているので、「2乗」というよりは、分散・共分散成分を知ることが大事であるから、
その期待値としては

$$
E[\frac{\partial L(w,\xi)}{\partial \xi_i} \frac{\partial L(w,\xi)}{\partial \xi_j}]
$$

となる。

これを
$$
g_{ij} = E[\frac{\partial L(w,\xi)}{\partial \xi_i} \frac{\partial L(w,\xi)}{\partial \xi_j}]
$$
と表す。



## フィッシャー情報量とリーマン計量の一致

前節で示したようにフィッシャー情報量は正方行列で

$$
g_{ij}=E[\frac{\partial L(w,\xi)}{\partial \xi_i} \frac{\partial L(w,\xi)}{\partial \xi_j}] = \int_\Omega \frac{\partial L(w,\xi)}{\partial \xi_i} \frac{\partial L(w,\xi)}{\partial \xi_j} p(w,\xi)d \mu
$$

一方、

$$
E[\frac{\partial}{\partial \xi_i} \frac{\partial}{\partial \xi_j} L(w,\xi)] = E[\frac{\partial}{\partial \xi_i}(\frac{1}{p(w,\xi)}\frac{\partial p(w,\xi)}{\partial \xi_j})]

$$
さらに変形して

$$
= E[-\frac{1}{p(w,\xi)^2}\frac{\partial p(w,\xi)}{\partial \xi_i})\frac{\partial p(w,\xi)}{\partial \xi_j})] + E[\frac{1}{p(w,\xi)}\frac{\partial}{\partial \xi_i}\frac{\partial p(w,\xi)}{\partial \xi_j}]
$$
第1項は
$$

- \int_\Omega (\frac{1}{p(w,\xi)} \frac{\partial p(w,\xi)}{\partial \xi_i}) (\frac{1}{p(w,\xi)} \frac{\partial p(w,\xi)}{\partial \xi_j}) d \mu
= - \int_\Omega \frac{\partial L(w,\xi)}{\partial \xi_i}\frac{\partial L(w,\xi)}{\partial \xi_j} d \mu
$$

第2項は
$$
\int_\Omega \frac{1}{p(w,\xi)}\frac{\partial}{\partial \xi_i}\frac{\partial p(w,\xi)}{\partial \xi_j} p(w,\xi) d \mu = \int_\Omega \frac{\partial}{\partial \xi_i}\frac{\partial p(w,\xi)}{\partial \xi_j} d\mu = 0
$$

なので、結局、

$$
E[\frac{\partial}{\partial \xi_i} \frac{\partial}{\partial \xi_j} L(w,\xi)] = - \int_\Omega \frac{\partial L(w,\xi)}{\partial \xi_i}\frac{\partial L(w,\xi)}{\partial \xi_j} d \mu
$$
つまり、フィッシャー情報量は、対数尤度関数の2階の微分に(-1)を掛けたものになっていることが解る。

2階の微分というのは、曲がり具合〜曲率であるから、この行列は多様体の曲率に相当する量であり、リーマン計量のことである。

ここで、推定のよさであるフィッシャー情報行列がリーマン計量であることがわかった。

### リーマン計量と「内積」

リーマン計量というのは、局所にベクトルがあったときに、それが作る内積をどのくらいの大きさにするかを決める行列である。

いわゆる内積を返すようなリーマン計量が単位行列であり、曲がっていれば、単位行列ではなくなることや、曲がっているところの「長さ」を測るときに用いられた行列$G$がこのリーマン計量である。


## $p(w,\xi)$ と $\log{p(w,\xi)}$ との一般化

関数とその対数とを「一般化」する必要などあるとも感じられないが、この2つを取ることの意味・意義について腑に落ちるためには、この「一般化」について知っておくことは有用である。

フィッシャー情報量を別の形で表現する、ということをとっかかりに、この「一般化」を眺めてみることにする。

### フィッシャー情報量と球面と自然な計量

フィッシャー情報量は以下のようにも表せる。

$$
g_{ij}= \int_\Omega \frac{\partial L(w,\xi)}{\partial \xi_i} \frac{\partial L(w,\xi)}{\partial \xi_j} p(w,\xi)d \mu = \int_\Omega \frac{1}{p(w,\xi)}\frac{\partial p(w,\xi)}{\partial \xi_i}\frac{\partial p(w,\xi)}{\partial \xi_j} d \mu
$$ 

この式では

$\farc{\partial L(w,\xi)}{\partial \xi_i} p(w,\xi) = \frac{1} {p(w,\xi)}\frac{\partial p(w,\xi)}{\partial \xi_i}} p(w,\xi)$ とすることで、分母分子に現れる$p(w,\xi)$をキャンセル・アウトできる。

今、二つの$\frac{\partial L(w,\xi)}{\partial \xi_i}$ $\frac{\partial L(w,\xi)}{\partial \xi_j}$を不平等に扱って$p(w,\xi)$をキャンセルアウトする代わりに、2つを平等に使ってキャンセル・アウトすることを考える。


そのように考えると、以下も成り立つことが示せる。


$$
g_{ij} = \int_\Omega \frac{\partial 2\sqrt{p(w,\xi)}}{\partial \xi_i}\frac{\partial 2\sqrt{p(w,\xi)}}{\partial \xi_j} d \mu= \int_\Omega \frac{1}{p(w,\xi)}\frac{\partial p(w,\xi)}{\partial \xi_i}\frac{\partial p(w,\xi)}{\partial \xi_j} d \mu
$$
 
これは、$2\sqrt{p(w,\xi)}$という変換をすると、半径2の球面になって、そこでの計量は一定になる、ということに相当する。

$\Omega$が有限個の状態であるときには、その個数次元の球面を考えればよい。

### 不均等なペアに分解する

均等に分解することも可能だが、$p(w,\xi)$をキャンセル・アウトするだけであれば、色々な不均等分解が可能である。

次のような不均等分解が可能である。

$$
f^{(k)}(w,\xi)= k p(w,\xi)^\frac{1}{k}
$$

とおけば

$$
g_{ij} = \int_\Omega \frac{\partial f^{(k)}(w,\xi)}{\partial \xi_i}\frac{\partial f^{(\frac{1}{k})}(w,\xi)}{\partial \xi_i} d \mu
$$

が、そのような不均等分解であることが確かめられる。

ちなみに、このような$k$を用いた関数ではなく

$$
f^{(\alpha)}(w,\xi)= \frac{2}{1-\alpha} p(w,\xi)^\frac{1-\alpha}{2}
$$
と表現するのが通例である。

このような$\alpha$での表現にすることで、

$\alpha=0$のときに

$$
f^{(\alpha=0)}(w,\xi) = 2 p(w,\xi)^{\frac{1}{2}}
$$

となり、これは、「均等分解」のときの表現となる。


そして、
$$
f^{\alpha=-1}(w,\xi) = p(w,\xi)
$$

なわけであるが、このときのペアはもともと$\log(p(w,\xi))$であったが、
$$
f^{(\alpha=1)}(w,\xi) =\frac{2}{1-1} p(w,\xi)^{\frac{1-1}{2}}
$$
となって、計算できない式になる。
逆に
$$
f^{(\alpha=1)} = \log{p(w,\xi)}
$$
と定義することで、首尾一貫する。

この$f^{(\alpha)}(w,\xi)$と$f^{(-\alpha)}(w,xi)$とは、お互いに補い合って、フィッシャー情報量(リーマン計量)を分解する関係であり、双対になっているという。

### 平坦と直交と$\alpha= \pm 1$

$f^{(\alpha = \pm 1)}(w,\xi)$は、ある特別なパラメタの取り方をすると、「平坦」になることが知られている。

そのパラメタの取り方が

$$
\log{p(w)} = C(w) + F \theta - \phi(\theta)
$$
$$
p(w) = H(w) \eta
$$
という取り方である。

すべての$p(w)$がこのようにパラメタを取れるわけではないが、非常に多くの、そして確率モデルとして大事な場合にこのようなパラメタを取ることができることは知られており、そのような分布を、指数型分布族と呼ぶ。

そして、ここでは$F,H$という二つの行列を使って表していることからわかるように、$\theta$,$\eta$はそれぞれ「平坦」であるけれども、それを実現するための行列$F,H$はお互いに関係がない。

このようなときには、$\theta$,$\eta$座標系は直交していない。ただ、相互に平坦なだけである。

$H=F^{-1}$とすると、平坦であり、かつ、直交した2つの座標系にすることができる。

ちなみに、$F,H$などを使って、線形変換しても、平坦なままであることを、これらの接続が「アフィン接続」である、と呼ぶ理由である。


アフィン変換は、いわゆる行列による一次線形な変換に平行移動を加えたものであるが、ここで考えている$\alpha$を用いた関数が定義する接続は、アフィン接続であることが解っている。

$\alpha=\pm 1$の場合というのは、その中で特殊なものであって、「うまくパラメタを取れば平坦になる」ものである。

そして、アフィン接続なので、うまく変換することで、相互に直交なパラメタの取り方も存在することから、いつでも、平坦かつ直交な$\theta$,$\eta$の座標系のペアが取れる。

そして平坦ではないが、$\alpha=0$の場合というのは、2つの双対な座標系を取る、というこの仕組みの中で、「自己双対」である、という特殊な場合であること。

# 結論2

結局、尤度関数として我々が普通に使うものは、指数型分布族なので、必ず$\theta$,$\eta$座標系がとれて、それらは両方とも平坦で相互に直交している。ということが解った。

その背景には、異なる座標系での直線があり、異なる座標系にある二つの直線の長さ(KLd)をピタゴラスの定理的に計算できること、片方の座標系の部分空間へ、もう片方の座標系の垂線を引くことが簡単になること、などがあることもわかった。

そして、また、そこには、ユークリッド空間では、2つの座標系が、自己双対なために、1つの座標系で出来ていると見えること

2016-11-27 統計多様体

[][][]統計多様体 絵にしてみる

  • 参考pdf
  • 統計多様体は、確率密度分布(確率質量分布)の集合である。集合の要素は確率密度分布に相当する。確率密度分布がパラメタで表されているとき、点はパラメタが決めた座標に応じて並んでいる。その並んだものが統計多様体
  • たとえば、指数分布というものがある。1パラメタで決まる確率密度分布である。
  • これを並べると、パラメタが正の実数を取るので、多様体は、正の実数を埋め尽くしている。
  • とは言え、それはただの、半直線。
  • ここで、個々の点が、確率密度分布であることを思い出すために、個々の点に「確率密度分布の姿」をつけて表示しよう。
  • 指数分布は0以上の値に対して確率密度が定まっているから、それを第二軸とすれば、
  • 第一軸に指数分布のパラメタ、第二軸に指数分布の台に相当する値、第三軸に確率密度を取らせて

f:id:ryamada22:20161127111304p:image

# 統計多様体

lambda <- seq(from=0,to=1,length=50)
R <- seq(from=0,to=3,length=50)
lambda.R <- as.matrix(expand.grid(lambda,R))
v <- matrix(0,length(lambda),length(R))

for(i in 1:length(lambda)){
	v[i,] <- dexp(R,lambda[i])
}

library(rgl)
plot3d(cbind(lambda.R,c(v)),xlab="parameter",ylab="x",zlab="prob")
  • 2パラメタ確率密度関数の場合に同様のことをしようとすると4次元必要だが、確率密度を色にすれば

f:id:ryamada22:20161127112656p:image

mu <- seq(from=-2,to=2,length=10)
sigma <- seq(from=1,to=2,length=10)
R <- seq(from=-3,to=3,length=10)

m.s.R <- expand.grid(mu,sigma,R)
v <- rep(0,length(m.s.R[,1]))
for(i in 1:length(m.s.R[,1])){
	v[i] <- dnorm(m.s.R[i,3],m.s.R[i,1],m.s.R[i,2])
}

plot3d(m.s.R,xlab="mean",ylab="s",zlab="x")
col <- v/max(v)
spheres3d(m.s.R,radius=0.1,color=rgb(col,1-col,1))

[][][]統計多様体 どうして「期待値」を出すの?

  • 統計多様体のフィッシャー情報行列を出しましょうというと
    • E¥[¥partial_{¥theta_i} ¥log(p_¥theta (w) ¥partial_{¥theta_j}¥log_¥theta (w) ¥]とか、『期待値』を計算する
    • どうしてかっていうと、『各点』は確率密度分布で、それぞれの『点』である確率密度分布と、その近傍の『点』である確率密度分布とをつないだときの変化具合は、「分布」と「分布」の違いで見るわけだが、それぞれの分布が台全体に広がっているので、「違い」をスカラーにするには、全体を見渡す必要があって、その見渡すときには、発生確率で重みづけた平均(期待値)が適当だから
  • じゃあ、どうして2階偏微分か、というと、「傾き」に興味があるのではなくて、「傾きがどれくらい曲がっているか・変化しているか」に興味があるから