2010-12-05

相関のある２つの確率変数の同時分布　さらに続き

R 分布確率密度分布

こちらの続き
多変量標準正規分布
- 複数の確率変数( $X_1,X_2,...,X_k$ )があり、それらは相互に独立
- 個々の確率変数は、平均０、分散１の正規分布
  - $X_i \sim N(0,1)$
  - $P(x_i)=\frac{1}{\sqrt{2\pi}} e^{-\frac{x_i^2}{2}}$
  - 原点からの距離の２乗 $|x|^2$ に応じて、確率が一定割合で減じる分布
- 複数のが独立であるときの同時確率分布
  - 同時確率は個々の変数の確率の積
  - $P(\bf{x}=(x_1,x_2,...,x_k)) = (\frac{1}{\sqrt{2\pi}})^k e^{-\frac{\sum_{i=1}^{k} (x_i^2)}{2}}$
  - $P(\bf{x})= (\frac{1}{\sqrt{2\pi}})^k e^{-\frac{|\bf{x}|^2}{2}}$
  - 距離の２乗が指数になっている関数(正規分布)は以下のことから、同時分布において、特徴的な挙動をする
    - 個々の確率変数の分布が距離の２乗に応じており、距離の２乗が指数関数の指数になっていること、独立変数の同時確率は、各変数の確率の積であること、積は指数の和であること、ユークリッド距離は $\bf{x}^2=\sum_{i=1}^k x_i^2$ であることから、同時確率にしても、「距離」( $|\bf{x}|^2$ )を指数とした関数形であることがわかる
- 多変量(k変量)標準正規分布
  - 同確率等高線はk次元球を描く( $|\bf{x}|$ が等しい点の集合が球面)
  - 変量間に相関があるときは、k次元楕円(楕球)を描く
- 多次元球と多次元楕球
  - k次元球面上の点は、k個の相互に直交する軸への射影の長さを座標として表せる(正規直交座標)
  - 正規直交座標軸は、k次元空間にて回転することができる
  - 多次元球を正規直交座標の各軸について拡縮したものが多次元楕球であり、その逆も言える
  - $\sum_{i=1}^k x_i^2=C$ が球
  - $\sum_{i=1}^k (\frac{x_i}{a_i})^2=C$ が楕球
  - 回転を表す行列は、回転後の正規直交座標軸を表す単位ベクトル $v_i$ を列ベクトルとしたもの( $(\bf{v_i})$ )
  - 軸方向への拡縮を表す行列 $L$ は、対角成分のみの行列で、対角成分は正。第i番対角成分の値は、第i番軸の拡縮率
  - は楕円座標、は球座標
    - $E_c=L V S_p$
    - $S_p=V^{-1} L^{-1}E_c$
- 相関のある多次元(k変量)標準正規分布
  - $\bf{x}$ という観察がなされる確率密度の等高線が楕球 $E_c$ で描かれているという
  - その確率密度 $Pr(\bf{x}|V,L)$ は $x'=V^{-1} L^{-1}x$ の $x'$ が相関のない多次元(k次元)標準正規分布でのそれ( $Pr(\bf{x'}|N_k)$ )である
  - $Pr(\bf{x'}|N_k)=\prod_{i=1}^k Pr(x_i')$ である
正規分布と任意分布とを対応させる
- $X$ は標準正規分布に従うものとする
- $Y$ は何かしらの分布に従うものとする
- それぞれの累積確率密度がなる関係にあるとする
  - $x=(Q_{norm})^{-1}(Q_{fx}(y))$ ということ
- それぞれの確率密度分布は $Pr_{norm}(x)=\frac{d Q_{norm}(x)}{dx}$ , $Pr_{fx}(y)=\frac{d Q_{fx}(y)}{dy}$ と表すとする
- $Pr_{norm}(x)=\frac{d Q_{norm}(x)}{dx} = \frac{d Q_{fx}(y)}{dx} = \frac{d Q_{fx}(y)}{dy} \frac{dy}{dx}=Pr_{fx}(y)\frac{dy}{dx}$
- したがって、 $\frac{dx}{dy} = \frac{Pr_{fx}(y)}{Pr_{norm}(x)}$
多変量(k変量)任意分布
- $X_i$ は標準正規分布に従うものとする
- $Y_i$ は何かしらの分布に従うものとする
- それぞれの累積確率密度が $Q_{norm}(\bf{x})=Q_{fx}(\bf{y})$ なる関係にあるとする
- それぞれの確率密度分布は $Pr_{norm}(\bf{x})=\frac{d Q_{norm}(\bf{x})}{d\bf{x}}$ , $Pr_{fx}(\bf{y})=\frac{d Q_{fx}(\bf{y})}{d\bf{y}}$ と表すとする
- $Pr_{norm}(\bf{x})=\frac{d Q_{norm}(\bf{x})}{d\bf{x}} = \frac{d Q_{fx}(\bf{y})}{d\bf{x}} = \frac{d Q_{fx}(\bf{y})}{d\bf{y}} \frac{d\bf{y}}{d\bf{x}}=Pr_{fx}(\bf{y})\frac{d\bf{y}}{d\bf{x}}$
- したがって、 $Pr_{fx}(\bf{y})=Pr_{norm}({\bf{x}}) \frac{d\bf{x}}{d\bf{y}}$
- $Pr_{fx}(\bf{y}|V,L)=Pr_{norm}({\bf{x}}|V,L) \frac{d\bf{x}}{d\bf{y}} =Pr_{norm}({\bf{x'}}|N_k) \frac{d\bf{x}}{d\bf{y}}$
- $Pr_{fx}(\bf{y}|V,L)=\prod_{i=1}^k Pr_{norm} x'_i \frac{dx_i}{dy_i}$
- $Pr_{fx}(\bf{y}|V,L)= \prod_{i=1}^k Pr_{norm} x'_i \frac{Pr_{fx}(y)}{Pr_{norm}(x)}$
- この等式は、少なくとも、 $\bf(x)=\bf(x')$ (独立なとき)には成立しているが、…この式変形でよいかな…
仮のRソース

multExp<-function(ys,lambdas,R){
	if(!is.matrix(ys))ys<-matrix(ys,nrow=1)
	# 特異値分解
	ER = sin(pi*R/2)
	ER.e=eigen(ER)
	if(min(ER.e[[1]])<0){
		print("Non-positive eigenvalue(s)!")
		print("Calculation was forced by non-zero eigenvalue(s) to 0.000000001.")
	}
	ER.e.ori=ER.e
	ER.e[[1]][ER.e[[1]]<0]=0.000000001
	# 標準正規分布の累積確率対応座標
	xs.el<-ys
	for(i in 1:length(ys[1,])){
		xs.el[,i]<-qnorm(pexp(ys[,i],1/lambdas[i]))
	}
	# 楕円座標の正球化
	# 固有値が0だと困る…
	# まず、楕球を回転させて、拡縮軸が標準化軸になるようにして (solve(ER.e[[2]])%*%t(xs.el))) 
	# 次いで、拡縮 (diag(1/sqrt(ER.e[[1]]))
	# 軸を元に戻す ER.e[[2]]
	xs.sp<-t(ER.e[[2]]%*%(diag(1/sqrt(ER.e[[1]]))%*%(solve(ER.e[[2]])%*%t(xs.el))))
	
	# 正球座標での確率
	#Pr.sp<-exp(apply(dnorm(xs.sp,log=TRUE),1,sum))
	Pr.sp<-dnorm(xs.sp,log=TRUE)
	Pr.el<-dnorm(xs.el,log=TRUE)
	Pr.y<-dexp(ys,log=TRUE)
	for(i in 1:length(ys[1,])){
		Pr.y[,i]<-dexp(ys[,i],1/lambdas[i])
	}

	Pr<-exp(apply(Pr.sp-Pr.el+Pr.y,1,sum))
	Pr.sp<-exp(apply(Pr.sp,1,sum))
	Pr.el<-exp(apply(Pr.el,1,sum))
	Pr.y<-exp(apply(Pr.y,1,sum))
	return(list(ys=ys,lambdas=lambdas,R=R,ER.e.ori=ER.e.ori,ER=ER,xs.el=xs.el,xs.sp=xs.sp,Pr=Pr,Pr.sp=Pr.sp,Pr.el=Pr.el,Pr.y=Pr.y))
}
plotMultExp<-function(multExpOut,ax=c(1,2),slice=1,theta=80,phi=80){
	# 算出した生起確率を３変数次元のアレイに納める
	PrArray<-array(multExpOut$Pr,c(length(tmpx),length(tmpx),length(tmpx)))
	# 第１変数の値をアレイに納める
	TArray1<-array(ys[,ax[1]],c(length(tmpx),length(tmpx),length(tmpx)))
	# 第２変数の値をアレイに納める
	TArray2<-array(ys[,ax[2]],c(length(tmpx),length(tmpx),length(tmpx)))

	#persp(tmpx,tmpx,PrArray[,,1])
	# プロットする
	persp(TArray1[,1,slice],TArray2[1,,slice],PrArray[,,slice],theta=theta,phi=phi,col="green",xlab=ax[1],ylab=ax[2])
}

# 変数間の相関強度の行列
nV<-3
R<-matrix(0,nV,nV)
R[1,2]<-R[2,1]<-0.4
R[1,3]<-R[3,1]<-0.3
R[2,3]<-R[3,2]<-0.2
diag(R)<-1
tmpx<-seq(from=0,to=1,length.out=20)
tmpx[1]<-0.001
tmpy<-tmpx
tmpz<-tmpx
#tmpz<-c(1)
#tmpx<-1:10
ys<-as.matrix(expand.grid(tmpx,tmpy,tmpz))
ysOrder<-as.matrix(expand.grid(1:length(tmpx),1:length(tmpy),1:length(tmpz)))

lambdas<-c(3,2,1)
multExpOut<-multExp(ys=ys,lambdas=lambdas,R=R)
multExpOut$ER.e
Pr<-multExpOut$Pr

Pr.Array<-Pr.sp.Array<-Pr.el.Array<-Pr.y.Array<-array(0,c(length(tmpx),length(tmpy),length(tmpz)))
for(i in 1:length(Pr)){
	Pr.Array[ysOrder[i,1],ysOrder[i,2],ysOrder[i,3]]<-Pr[i]
	Pr.sp.Array[ysOrder[i,1],ysOrder[i,2],ysOrder[i,3]]<-multExpOut$Pr.sp[i]
	Pr.el.Array[ysOrder[i,1],ysOrder[i,2],ysOrder[i,3]]<-multExpOut$Pr.el[i]
	Pr.y.Array[ysOrder[i,1],ysOrder[i,2],ysOrder[i,3]]<-multExpOut$Pr.y[i]
}
logPr.Array<-log(Pr.Array)

zlim<-c(min(logPr.Array),max(logPr.Array))
for(i in 1:length(tmpz)){
	tmpi<-100+i
	file.name=paste("test6",tmpi,".jpg",sep="")
	jpeg(file=file.name)
	filled.contour(logPr.Array[,,i],zlim=zlim)
	#persp(PrArray[,,i],col="green",phi=60,theta=70)

	dev.off()
}


library(rgl)
open3d()
logPr<-log(Pr)
plot3d(ys[,1],ys[,2],ys[,3],type="s",col=gray((max(logPr)-logPr)/(max(logPr)-min(logPr))))