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

相関のある２つの確率変数の和の分布

R 分布 確率密度分布
  • かなり怪しい話。
  • こちらとこちらで、独立する確率変数の和の分布について書いた
  • 和の分布の確率密度分布は、２つの変数の分布の積であって、和の分布の特性関数は、２つの変数の分布の特性関数の積であること、何かしら、うまい工夫をすると、相関する２確率変数の和の分布に関して、解析的に処理できる部分がないのか、それには、２つの複素変数の(２つの複素平面)の間に直交でない関係を入れればよいのではないか、というようなことを、とりとめなく書いた
  • 特性関数まで持っていくのは大変なので、まずは、確率密度関数のレベルでどうなるかを考える
  • 平面に原点を同じくする２つの座標系(XYとZW)を置く
  • ZW座標系はいわゆるデカルト座標とする。その座標を(z,w)で表す
  • XY座標系はZ軸を中心に対称なX軸とY軸とからなる
  • X軸はZ軸から-\theta回転したもの、Y軸は\theta回転したものとし、X、Yは相互に2\thetaの角をなしている。\theta=\frac{\pi}{4}のとき、X,Yは直交している
  • このとき、ある点をXY座標系、ZW座標系で表した座標(x,y),(z,w)とは、次の行列によって、変換される
  • \begin{pmatrix} z \\ w \end{pmatrix} = \begin{pmatrix} \cos(\theta) & \cos(\theta) \\ -\sin(\theta) & \sin(\theta) \end{pmatrix} \begin{pmatrix} x \\ y \end{pmatrix}
  • \begin{pmatrix} x \\ y \end{pmatrix} = \begin{pmatrix} \cos(\theta) & \cos(\theta) \\ -\sin(\theta) & \sin(\theta) \end{pmatrix}^{-1} \begin{pmatrix} z \\ w \end{pmatrix} = \begin{pmatrix} \frac{1}{2\cos(\theta)} & \frac{-1}{2\sin(\theta)} \\ \frac{1}{2\cos(\theta)} & \frac{1}{2\sin(\theta)} \end{pmatrix}\begin{pmatrix} z \\ w \end{pmatrix}
  • この座標変換をRでプロットするのが次のソース
par(mfcol=c(2,2))
ttt<-seq(from=0,to=2*pi,length.out=100)
xxx<-cos(ttt)
yyy<-sin(ttt)
xlim<-ylim<-c(-2,2)
k<-0.3 # 相関の強さ k=1は無相関、k=0は完全一致
t<-pi/4 *k  　
Mxy2ab<-matrix(c(cos(t),cos(t),-sin(t),sin(t)),2,2,byrow=TRUE)
zzz<-www<-xxx
for(i in 1:length(xxx)){
 tmp<-Mxy2ab%*% matrix(c(xxx[i],yyy[i]),2,1)
 zzz[i]<-tmp[1]
 www[i]<-tmp[2]
}
plot(xxx,yyy,xlim=xlim,ylim=ylim)
plot(zzz,www,xlim=xlim,ylim=ylim)
minmin<--10
maxmax<-10
ttt<-seq(from=minmin,to=maxmax,length.out=100)
xxx<-ttt
yyy<-xxx
xlim<-ylim<-c(2*minmin,2*maxmax)
zzz<-www<-xxx
for(i in 1:length(xxx)){
 tmp<-Mxy2ab%*% matrix(c(xxx[i],yyy[i]),2,1)
 zzz[i]<-tmp[1]
 www[i]<-tmp[2]
}
plot(xxx,yyy,xlim=xlim,ylim=ylim)
plot(zzz,www,xlim=xlim,ylim=ylim)
par(mfcol=c(1,1))
  • ２つの確率変数X'、Y'が独立である場合をまず考える
    • X'、Y'がそれぞれ、値x',y'をとるときその確率はPr(x',y')=f(x')\times g(y')で与えられる
    • Pr(x',y')=C(Cは定数)を満足する点が平面上に等高線を描く
    • たとえば、２つの確率変数同じであれば、
  • 点P=(x,y)=(z,w)は、確率変数Xが値xをとり、かつ、確率変数Yが値yを取る事象に対応するものとする
    • 点Pの
  • 以下はそのためのRソースメモ
  • 掲載図は、正規分布と指数分布とが、ある程度の相関があるときの確率密度分布。緑ののっぺりした領域は指数分布が非負領域のみの分布であるために確率が０である領域。
  • 掲載図の横軸は、２つの確率変数の和に比例した軸なので、和がある値であるとき確率は、垂直軸に関して、積分した値となる、そんな風に描出してある
minXY<--10
maxXY<-10
kizami<-0.1
xlim<-ylim<-c(-2*minXY,2*maxXY)
x<-seq(from=-maxXY,to=maxXY,by=kizami)
y<-x
XY<-as.matrix(expand.grid(x,y))

k<-0.3 # 相関の強さ k=1は無相関、k=0は完全一致

t<-pi/4 *k  　
Mxy2ab<-matrix(c(cos(t),cos(t),-sin(t),sin(t)),2,2,byrow=TRUE)

AB<-matrix(0,length(XY[,1]),2)
for(i in 1:length(AB[,1])){
 tmp<-Mxy2ab %*% matrix(c(XY[i,]),2,1)
 AB[i,]<-tmp
}

plot(AB[,1],AB[,2],xlim=xlim,ylim=ylim)

# z = x + y = a /cos(t)

Z<-AB[,1]/cos(t)

library(rgl)
plot3d(AB[,1],AB[,2],Z)

##############

A<-B<-seq(from=-maxXY,to=maxXY,by=kizami)
Z<-A/cos(t)

Zmat<-matrix(rep(Z,length(B)),length(Z),length(B),byrow=FALSE)


persp(Zmat,col="#789ABC")

Mab2xy<-matrix(c(1/(2*cos(t)),-1/(2*sin(t)),1/(2*cos(t)),1/(2*sin(t))),2,2,byrow=TRUE)

X<-Y<-Zmat
for(i in 1:length(Zmat[,1])){
	for(j in 1:length(Zmat[1,])){
		tmp<-Mab2xy%*% matrix(c(A[i],B[j]),2,1)
		X[i,j]<-tmp[1]
		Y[i,j]<-tmp[2]
	}
}

#Xpr<-dnorm(X,mean=0,sd=1)
Xpr<-dexp(X,0.01)
Ypr<-dnorm(Y,mean=0,sd=5)
#Ypr<-Xpr
Zpr<-Xpr*Ypr
par(mfcol=c(2,2))
image(Xpr)
image(Ypr)
image(Zpr)


#persp(Zpr,col="#789ABC")
#image(Zpr,col=)
filled.contour(Zpr, color = terrain.colors, asp = 1)

plot3d(row(Zpr),col(Zpr),Zpr,type="h",col=rainbow(100))