2017-03-05

双対平坦座標系が双直交する

具体例がないとわかりにくいので、具体例から入る
正規分布 $p= \frac{1}{\sqrt{2\pi s^2}}e^{-\frac{(x-m)^2}{2s^2}}$ を $\theta$ 座標系と $\eta$ 座標系で表すことができる
$\theta$ 座標系というのは、分布関数を指数族表現したときに現れるパラメタであって、正規分布の場合には $\log{p} = \theta_1 x + \theta_2 x^2 + ...$ となるようなそれである
具体的には $\theta_1 = \frac{m}{s^2},\theta_2 = -\frac{1}{2s^2}$
他方、 $\eta$ 座標系というのは、指数族表現して現れた $\theta$ パラメタに掛け合わされる $x$ の関数の期待値を取ることになっているから、 $x,x^2$ の期待値のこと
正規分布の場合、 $x$ の期待値は、 $m$ だし、 $x^2$ の期待値は $(x-m)^2=s^2$ なる関係があるから、 $m^2+s^2$ が期待値である。 $\eta_1=m,\eta_2=m^2+s^2$ というパラメタを取る、と言うことである
$\theta$ 座標系、 $\eta$ 座標系のそれぞれで、接ベクトル(関数を偏微分したもの)を取り、その関係を調べると「双直交」だとわかるわけであるが
$\theta$ 、 $\eta$ 座標系ではそれぞれ $\log{p}$ , $p$ を偏微分することになっている
それは、フィッシャー情報量が $g_{ij} =\int \frac{\partial \log{p}}{d\partial \xi_i}2\sqrt{p}}{\partial \xi_j} dx$ と表されることと関係していて、もっと一般的には $g_{ij} = \int \frac{\partial l^{(\alpha)}}{\partial \xi_i}\frac{\partial l^{(-\alpha)}}{\partial \xi_j}dx$ ( $l^{(\alpha)} = \frac{2}{1-\alpha} p^{\frac{1-\alpha}{2}}$ と表されるときの $\alpha=1$ が[l^{(1)}=\log{p}]であって、そのときに現れるパラメタが $\theta$ 、 $\alpha=-1$ が $l^{(-1)}$ であって、そのときに現れるパラメタが $\eta$ であることと結びついている
別の説明もできる。 $\theta,\eta$ のどちらにも $\log{p}$ の偏微分を取るのだが、偏微分同士の積の期待値を取ろうとすると、 $\int \frac{\partial \log{p}}{\theta_i}\frac{\partial \log{p}}{\eta_j}\times p dx$ となるが、 $\frac{\partial \log{p}}{\eta_j} = \frac{1}{p}\frac{\partial p}{\eta_j}$ であるから、 $\times p$ がキャンセルアウトする、とも見える
実際、この二つの座標系について、それぞれのパラメタを微小変化させ、それに対応する正規分布(とその $\log{p}$ )との微小変化を計算し、偏微分相当の値を算出し、それを無限次元ベクトル(の有限離散版)として、 $\frac{\partial \log{p}}{\partial \theta_i}$ と $\frac{\partial p}{\partial \eta_j}$ との内積を取ると、以下に示すように、単位行列が現れ、『双直交』であることの例示ができる

# 座標変換関数
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.thata2eta <- 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])
}

# 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
dx <- 1/1000
x <- seq(from=-10,to=10,by=dx)
p <- dnorm(x,m0,s0)
plot(x,p)
# 対応する$\theta,\eta$を求める

thetas <- my.ms2theta(m0,s0)
etas <- my.ms2eta(m0,s0)
# $theta_1$の偏微分を求める。その増減の様子を描く

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$についても同様に求める

d.theta2 <- 0.0001

d.eta1 <- 0.0001
d.eta2 <- 0.0001
dp.dtheta1 <- (my.dnorm.theta(x,thetas[1]+d.theta1,thetas[2]) - my.dnorm.theta(x,thetas[1],thetas[2]))/d.theta1

dp.dtheta2 <- (my.dnorm.theta(x,thetas[1],thetas[2]+d.theta2) - my.dnorm.theta(x,thetas[1],thetas[2]))/d.theta2

dp.deta1 <- (my.dnorm.eta(x,etas[1]+d.eta1,etas[2]) - my.dnorm.eta(x,etas[1],etas[2]))/d.eta1
dp.deta2 <- (my.dnorm.eta(x,etas[1],etas[2]+d.eta2) - my.dnorm.eta(x,etas[1],etas[2]))/d.eta2
# 双直交性の確認 (1,0,0,1)となる

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]