2012-06-30

確率の比

多次元空間に２つの多次元(標準)正規分布があるとする
この空間の点において、２つの正規分布の確率p1,p2の比を求めることにする
もちろん２つの正規分布の中心からの距離に応じて、それぞれの確率を計算してからその比を取ってよい
何かしらの理由でそうできないとする
多次元正規分布の中心から距離(距離の二乗)でまとめたときの分布はカイ分布(カイ二乗分布)
多次元正規分布の中心でない点からの距離(距離の二乗)でまとめたときの分布は非心カイ分布(非心カイ二乗分布)
非心カイ二乗分布の比がどうなるかを考えてみる
こちらにあるように、非心カイ二乗分布はという
- ここで $I_{\frac{k}{2}-1}(\sqrt{\lambda x})$ はmodified bessel function of the first kind
別の書き方もある
- ただし $_0F_1(;\frac{k}{2},\frac{\lambda x}{4})$ は超幾何関数、 $P(k,x)$ は自由度kのカイ二乗分布の確率密度関数
２つの非心カイ二乗分布の比を取るときには、が２つ異なるときに、の点で比較する
- このとき $I_{\frac{k}{2}-1}(\sqrt{\lambda x})$ は両分布でともに $I_{\frac{k}{2}-1}(0)$ と同じになるので、、比を取ると消えるかと思えば、 $x=0$ のときは $\frac{0}{0}$ となるので、そうはいかない
- 他方、超幾何関数の方は、ともに $x=0$ でともに0でない値に収束し、その値が異なる $\lambda$ で等しいので、超幾何関数を用いて比を取るときには、打ち消しあう
- $x^t$ の項はどちらでも消える
やってみる(多次元正規分布の確率の比とも比較する)
２つのk次元正規分布があるときに、ある点での２つの分布の生起確率の比は、その点における、２つの多次元正規分布に関して得られる自由度kの非心カイ自乗分布の比のカイ自乗値0における値(自由度が３以上では $\frac{0}{0}$ になるので、極限)と同じ(らしい)

# k次元。中心はxy平面にある。xy平面上の点だけを見る
k <-4
c1 <- rep(0,k)
c2 <- c(1,rep(0,k-1))
x <- y <- seq(from=-2,to=3,by=pi/30)
xy <- expand.grid(x,y)
r1 <- apply(xy^2,1,sum)
r2 <- apply(matrix(c(xy[,1]-1,xy[,2]),ncol=2)^2,1,sum)
p <- exp(-(r1-r2)/2) 
m.p <- matrix(p,ncol=length(x))
image(m.p)
contour(m.p)
library(rgl)
plot3d(as.matrix(cbind(xy,p)))
q <- exp(-1/2*(r1-r2))
plot(p/q)
m.q <- matrix(q,ncol=length(x))
image(m.q)
contour(m.q)
library(rgl)
plot3d(as.matrix(cbind(xy,q)))

ベッセル関数・超幾何関数の比の収束についてのメモが以下のソース

library(gsl)

k <- 3
x <- seq(from=0,to=10,length=10000)

lambda1 <- 2
lambda2 <- 3

b1 <- besselI(sqrt(lambda1*x),k/2-1)
b2 <- besselI(sqrt(lambda2*x),k/2-1)

b <- data.frame(x = x, b1=b1,b2=b2)
library(ggplot2)
group <- c(rep(1,length(b1)),rep(2,length(b2)))
gp <- ggplot(data = b, aes(x=rep(x,2),y=c(b1,b2),colour=factor(group)))

gp <- gp + geom_line()

gp

gp.2 <- ggplot(data = b, aes(x=x,y=c(b1/b2)))

gp.2 <- gp.2 + geom_line()

gp.2


h1 <- hyperg_0F1(k/2,lambda1*x/4)
h2 <- hyperg_0F1(k/2,lambda2*x/4)

h <- data.frame(x = x, b1=h1,b2=h2)
library(ggplot2)
group <- c(rep(1,length(b1)),rep(2,length(b2)))
gq <- ggplot(data = h, aes(x=rep(x,2),y=c(b1,b2),colour=factor(group)))

gq <- gq + geom_line()

gq

gq.2 <- ggplot(data = h, aes(x=x,y=c(b1/b2)))

gq.2 <- gq.2 + geom_line()

gq.2