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

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

2018-01-24 Maximum Mean DIscrepancy その3

[][][][]Maximum Mean DIscrepancy その3

  • Rのkernlab
library(kernlab)
# create data
x <- matrix(runif(300),100)
y <- matrix(runif(300)+1,100)


mmdo <- kmmd(x, y)

mmdo
> mmdo
Kernel Maximum Mean Discrepancy object of class "kmmd" 
 
Gaussian Radial Basis kernel function. 
 Hyperparameter : sigma =  1.18155987898832 


 H0 Hypothesis rejected :  TRUE
 Rademacher bound :  0.546163676520458

 1st and 3rd order MMD Statistics :  1.03692666322171 1.06709240109231


> mmdo@mmdstats
[1] 1.036927 1.067092
  • MMD統計量は2つ出る。mmdo@mmdstats
  • こちらのコードをたどると
mmd1 <- sqrt(max(0,sumKxx/(m*m) + sumKyy/(n*n) - 2/m/n* sumKxy))
mmd3 <- sum(hu)/M/(M-1)
    • とあり、Biased MMDと、Unbiased MMD(のN=M(標本数が同じ)の場合の計算式)になっている
  • Python

2018-01-23 フーリエ、巻きつける

[][][][]フーリエ、巻きつける

  • 分布関数があったときに、E¥[e^{tX}¥]としてやるとモーメント母関数。tXが実数tに関して無限の広がりを持つので、無限に複雑な分布関数の情報を担わせることができる
  • ちょっと変える。E¥[e^{itX}¥]。これは特性関数。複素数を導入することで、無限遠まで伸びる実数直線を単位円周に「巻きつけ」た。円周は有限だが、無限に細かく評価できるから、無限の情報を担わせることができる。
  • 二次元で考える。二次元平面のすべての方向について「巻きつける」と単位球が現れる。これに二次元平面の分布の情報は巻きつくだろうか…

2018-01-22 ぱらぱらめくる『カーネル法入門』

[][]ぱらぱらめくる『カーネル法入門』

  • 章末にまとめが書いてある。それをかいつまむと、この本の主張が構成に即してわかるように思うのでそのようにする
  • 1,2章
    • データの高次モーメントを取り出すことができるカーネル法
    • 特徴空間に取り出してやり、そこで線形解析手法を実行することで非線形解析が可能になる
    • どうしてそんなことが出来るかといえば、線形解析手法は内積計算の積み重ねであり、特徴空間では内積計算が、カーネルによって定まっているから
  • 3章
    • カーネルと特徴空間を使って、様々な線形解析手法の非線形版がある
  • 4章
    • カーネル法の代表であるサポートベクターマシンをその拡張に焦点を当てて扱う
  • 5章
    • カーネル法に代表される統計的学習を誤差の観点から扱う
  • 6章
    • 数学的側面
    • 正定値、負定値、フーリエ変換
  • 7章
    • カーネル法はカーネルさえ定義できれば、元のデータはユークリッド空間になくてもよい
    • ストリング構造(文字列)、グラフ
  • 8章
    • カーネル法にてデータ解析をするのは、確率変数を写像空間に持ち上げてそこで線形解析をすることだった
    • もっと根本に戻って、「確率変数」がどうなっているのかを特徴空間で考える。平均やモーメントが出てくる
    • 特性関数の一般化が得られる
    • 分布が写像空間で扱えるので、2分布の異同なども、検定等計量(U統計量)に落とせる
  • 9章
    • 分布を特徴写像で考えることの延長として、データの独立性、条件付独立性などの議論が行える。「分布」の台を一般化して考えることができる
  • 10章
    • 関数データ解析、スプライン平滑化、確率過程もカーネル法の枠組みで考えることができる

2018-01-21 特性関数・指数型分布族・情報幾何

[][][][][]特性関数・指数型分布族・情報幾何

  • メモ的なRmdファイル
---
title: "特性関数・指数型分布族・情報幾何"
author: "ryamada"
date: "2018年1月18日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

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

# 指数型分布族

```{r}
my.make.exponential <- function(C,F,A){
  list(C=C,F=F,A=A)
}
my.instance.exponential <- function(myexp,theta){
  list(C=myexp$C,F=myexp$F,A=myexp$A,theta=theta)
}
my.value.exponential <- function(myins,x){
  Cx <- myins$C(x)
  Fx <- matrix(0,length(x),length(myins$F))
  for(i in 1:length(Fx[1,])){
    Fx[,i] <- myins$F[[i]](x)*myins$theta[i]
  }
  Atheta <- myins$A(myins$theta)
  tmp <- Cx + apply(Fx,1,sum) - rep(Atheta,length(x))
  return(exp(tmp))
}
```

# 正規分布

```{r}
C.normal <- function(x){
  return(rep(log(1/sqrt(2*pi)),length(x)))
}
F.normal <- list()
F.normal[[1]] <- function(x){
  return(x)
}
F.normal[[2]] <- function(x){
  return(x^2)
}
A.normal <- function(theta){
  return(-theta[1]^2/(4*theta[2]) - 1/2 * log(-2 * theta[2]))
}

NRM <- my.make.exponential(C.normal,F.normal,A.normal)
mu <- 3
s <- 2
theta <- c(mu/s^2,-1/(2*s^2))
NRM.1 <- my.instance.exponential(NRM,theta)

x <- seq(from=-3,to=3,length=100)
(NRM.1$F[[1]])
```

```{r}
NRM.1.x <- my.value.exponential(NRM.1,x)
NRM.2.x <- dnorm(x,mu,s)
plot(NRM.1.x,NRM.2.x)
plot(x,NRM.1.x)
```

# モーメント母関数と特性関数とその元となる関数

モーメントは関数は $E\[e^{tX}\]$。

特性関数は$E\[e^{itx}\]$。

$e^{tx}$,$e^{itx}$は、確率変数。
$X$の台$x$に関する関数と見れば、$t,x$の2つの異なる変数が定める関数。

その期待値を取ると、$x$については集約され、$t$のみの関数となる。

## $e^{tx}$と$e^{itx}$

$e^{tx}$,$e^{itx}$は、どちらも、確率変数$X$から作られる次元が増えた関数。

$e^{tx}$は次元が1つ増える。
$e^{itx}$は次元が2つ(実と虚)増える。
ただし、$e^{itx}$のとる値はMod=1に限定しているので、実は、複素数の「角度」の次元分(1)が増えているだけ。違いは、この角度の次元は周期性があること。

結局、$e^{tx}$は実数無限空間という次元を1つ加えるのに対し、$e^{itx}$は単位円周に角度に相当する次元を1つ加えている。

```{r}
library(rgl)

x <- seq(from=-3,to=1,length=100)
t <- seq(from=-3,to=1,length=100)
m <- 2
s <- 3

Px <- my.value.exponential(NRM.1,x)
plot(x,Px,type="l")

tx <- expand.grid(x,t)
etx <- exp(tx[,1] * tx[,2])
eitx <- exp(1i * tx[,1] * tx[,2])

Mod(eitx)

plot3d(tx[,1],tx[,2],etx)

plot3d(tx[,1],tx[,2],Re(eitx))
plot3d(tx[,1],tx[,2],Im(eitx))

plot(Im(eitx),etx)

plot3d(tx[,1],tx[,2],Arg(eitx))


```

## $E\[e^{tx}\]$と$E\[e^{itx}\]$


```{r}
my.eitx <- function(t,nrm,x=seq(from=-100,to=100,length=10^5),I=TRUE){
  p <- my.value.exponential(nrm,x)
  ret <- rep(0,length(t))
  q <- 1
  if(I){
    q <- 1i
  }
  for(i in 1:length(ret)){
    ret[i] <- sum(exp(q * t[i]*x) * p)/sum(p)
  }
  return(ret)
}

```

```{r}
t <- seq(from=-2,to=2,length=100)
Eetx <- my.eitx(t,NRM.1,I=FALSE)
Eeitx <- my.eitx(t,NRM.1,I=TRUE)
```
```{r}
plot(t,Eetx)
plot3d(t,Re(Eeitx),Im(Eeitx))
```

# 指数型分布族のlog-potential 関数とキュムラント

パラメタ[tex:\theta]はベクトル、それを変化させるパラメタ[tex:t]もベクトル。

その変化ベクトル[tex:t]に対して十分統計量[tex:F(x)]の線形和変数に関して、log partition関数が、ただそれだけで、特性関数を決める。

ここで十分統計量の線形和変数というのは、[tex:\theta]が定めるある特定の指数型分布族確率分布に対応する。

結局、指数型分布族が作っている分布の多様体のそれぞれの上に特性関数が定まっており、特性関数は、十分統計量が決める任意の確率変数のモーメント母情報をすべて持っているので、確率分布の変化具合ももれなく決まる。

とくに、個々の十分統計量を分離して、その1次モーメントの期待値を考えることもできるが、それは、情報幾何における、[tex:\eta]座標に相当する。

逆に言うと、十分統計量の取り方が決まっていて、それの期待値を決めてやると([tex:\eta]座標を決めてやると)、分布自体が確定するし(この指数型分布族は座標系の自由度しか持たない)し、それを決めているのもlog partition関数という、[tex:\theta]によってのみ決まる関数であることもわかる。


$$
E[e^{i\sum t_j F_j(x)}|\theta] = \int e^{i \sum t_j  F_j(x)} exp^{C(x) + \sum F(x)_j \theta_j - \psi(\theta)}dx = e^{\psi(\theta + i t) - \psi(\theta)}
$$

```{r}
my.exp.cumulant <- function(t.mat,myins,I=TRUE){
  ret <- rep(0,length(t.mat[,1]))
  q <- 1
  if(I){
    q <- 1i
  }
  Atheta2 <- myins$A(myins$theta)
  for(i in 1:length(ret)){
    Atheta1 <- myins$A(myins$theta+q*t.mat[i,]+0*1i)
    ret[i] <- Atheta1 - Atheta2
  }
  return(exp(ret))
}
```
```{r}
NRM.1$A(NRM.1$theta+0*1i+1)
```
```{r}
t <- seq(from=-2,to=2,length=100)
t.mat <- cbind(t,rep(0,length(t)))
# ここで、t.matのうち、第一十分統計量だけを変化させ、第二十分統計量を固定してやれば、第一十分統計量がXそのものなので、確率変数Xの特性関数がlog-partition関数から出ることも示せる
Eetx2 <- my.exp.cumulant(t.mat,NRM.1,I=FALSE)
Eeitx2 <- my.exp.cumulant(t.mat,NRM.1,I=TRUE)
```
```{r}
plot(t,Eetx2)
```
```{r}
plot(Re(Eeitx),Re(Eeitx2))
```

2018-01-20 カーネル法?KL divergence?

[][]カーネル法?KL divergence?

  • メモ的なRmdファイル
---
title: "カーネル法?KL divergence?"
author: "ryamada"
date: "2018年1月17日"
output: 
  html_document:
    toc: true
    toc_depth: 6
    number_section: true
---

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

# 分布の違いを定量する

分布の違いを定量する方法として、カルバック-ライブラー ダイバージェンス(KLd)が知られている。

また、カーネル法に基づく定量法も知られている。

これらの関係を考えてみることにする。

# KLd と JSd

確率変数 $X,Y$ が台 $\Omega$に定められており、その確率密度関数が$P_X,P_Y$であるとする。

このとき、
$$
KLd(X|Y) = \int_\Omega \log{\frac{P_X}{P_Y}} d P_X
$$
## 対称性・非対称性

KLdは非対称である。
KLdを基礎として対称化したものにJSdがある。

$$
JSd(X,Y) = JSd(Y,X) = \frac{1}{2}(KLd(X|Y)+KLd(Y|X)) = \frac{1}{2} (\int_\Omega \log{\frac{P_X}{P_Y}} d P_X + \int_\Omega \log{\frac{P_Y}{P_X}} d P_Y)
$$
整理して書き直すと

$$
JSd(X,Y) = \frac{1}{2} \int_\Omega (\log{P_X}-\log{P_Y}) d(P_X - P_Y) = \frac{1}{2} \int_\Omega (\log{P_X}-\log{P_Y}) (P_X - P_Y) dw
$$
# カーネル平均に基づく定量

$X$ と$Y$とからのランダム標本が$Z_X,Z_Y$とが得られているときに

$$
<f_X,f_Y>_H = \frac{1}{|Z_X||Z_Y|}\sum_{u \in Z_X}^{|Z_X|} \sum_{v \in Z_Y}^{|Z_Y|} K(u,v)
$$
を持って、$X$と$Y$との違いの定量値とすることもある。

これは、$X,Y$をカーネル$K(u,v)$が定めるReproducing Kernel Hilbert Space関数($f_X=E[K(.,X)],f_Y=E[K(.,Y]$)に対応付け、その2つの関数のヒルベルト空間($H$)での内積である。

ただし、$K(u,v)$はカーネル関数であって、$K(u,v)=K(v,u)$であり、正定値性を持つなど、制約はあるが、おおざっぱに言うと「近ければ大きく、遠ければ小さい値を返す対称性関数」のことと思っておけばよい。

良く用いられるものにガウシアン・カーネルがあり、
$$
K(u,v) = e^{-\gamma |u-v|^2}; \gamma > 0
$$
と表される。

## 分布関数既知の場合

上記は、標本から定量しているが、分布関数が知られているときには以下のようになる。

$$
<f_X,f_Y>_H = \int_{\Omega_X} \int_{\Omega_Y} K(P_X,P_Y) d P_X d P_Y = \int_{\Omega_X} \int_{\Omega_Y} K(P_X,P_Y) P_X P_Y dw_x dw_y
$$

# JSdとカーネル平均内積との比較

JSdの定義式は分布関数が知られているものとして定義されているので、両者を比較するには、以下の2式を比較することが適切である。

$$
JSd(X,Y) =  \frac{1}{2} \int_\Omega (\log{P_X}-\log{P_Y}) (P_X - P_Y) dw
$$

$$
<f_X,f_Y>_H = \int_{\Omega_X} \int_{\Omega_Y} K(P_X,P_Y) P_X P_Y dw_x dw_y
$$

両者の大きな違いは、JSdが1重積分、カーネル平均内積が2重積分であることである。

## 2重積分を1重積分にする

カーネルが
$$
K'(x,y) = 0; if x \ne y
$$
であるときには、$X,Y$に関する2重積分は1重積分となる。

$$
<f_X,f_Y>_H = \int_{\Omega_X} K'(P_X,P_Y) P_X P_Y dw_x
$$

したがって、JSdというのは、カーネル関数が$K'(x,y) = 1 (x=y), K'(x,y)= 0 (x \ne 0)$
であるとき

$(\log{P_X}-\log{P_Y}) (P_X - P_Y)$を$\Omega$に関して積分するのがJSdであり、
$P_X P_Y$ を積分するのがカーネル平均内積である、という違いであることがわかる。

# 平滑化・分布推定の道具としてのカーネル関数

ガウシアン・カーネルは、分布推定の用に用いることがある。

ある地点の密度を推定するにあたり、近傍の情報を使う方法に相当する。

$\gamma$が小さければ小さいほど、より広い範囲の情報を使うことに相当し、$\gamma$が大きければ大きいほど狭い範囲の情報のみを使うことに相当する。

また、有限個の標本情報に基づいて分布推定する場合には、$\gamma$が小さければ小さいほど、推定分布はなめらかになり、大きければ大きいほど推定分布は複雑になる。


推定においては、標本数が多ければ多いほど、$\gamma$を大きくすることが適切になる。

逆に言うと、KLd, JSdの場合には、分布関数が既知であるときに定義されているが、標本が得られているときには、その標本にカーネルを適用するなどして、分布を推定し、その推定分布についてKLd,JSdを計算するという手順を踏むことになる点で、KLd,JSdもカーネルと無縁ではない。