2010-12-06

相関のある複数の確率変数の同時分布　一般的に

R 分布確率密度分布

こちらの続き
を共有する標準正規分布に従う個の確率変数の同時分布を考える
- $\bf{Y}=(Y_i)$ :標準正規分布に従う確率変数
- $y_i \in Y_i$
- $E_i$ : $Y_i$ が従う確率密度分布(いずれも標準正規分布)
- :の確率(累積)分布
  - $Q_i(t_i)=\int_{-\infty}^{t_i} D_{NL}(y_i) dy_i$
- はの相関係数行列である
  - $r_{ij}=Cor(Y_i,Y_j)$
- なるがあり
  - $\bf{M}=\bf{VL^2V^{-1}$ と特異値分解する
  - このとき $\bf{y'}=\bf{VL^{-1}V^{-1y}$ とする
- $\bf{E}(\bf{y}=(y_i)|\bf{R})$ : $\bf{R}$ の下での $\bf{Y}$ の同時確率分布
- - ただし $\bf{I}$ は単位行列であり $\bf{R}=\bf{I}$ は独立を意味する
任意の分布に従う、複数の確率変数の同時分布を考える
- $k$ :確率変数の数
- $\bf{X}=(X_i);i=1,2,...,k$ :確率変数
- $x_i \in X_i$
- :が従う確率密度分布
  - $\int_{-\infty}^{\infty} D_i(x_i) dx_i = 1$
- :の確率(累積)分布
  - $P_i(s_i) = \int_{-\infty}^{s_i} D_i(x_i) dx_i$
- :の相関係数行列
  - $r_{ij} = Cor(X_i,X_j)$
- $\bf{D}(\bf{x}=(x_i)|\bf{R})$ : $\bf{R}$ の下での $\bf{X}$ の同時確率分布
ととの関係
- - $y_i=Q_i^{-1}(P_i(x_i))$
ととの関係を用いてを表す
- $\bf{D}(\bf{x}|\bf{R})=\bf{E}(\bf{y}|\bf{R})\prod_{i=1}^k \frac{D_i(x_i)}{E_i(y_i)}$
- $\bf{D}(\bf{x}|\bf{R})=\bf{E}(\bf{y'}|\bf{R}=\bf{I})\prod_{i=1}^k \frac{D_i(x_i)}{E_i(y_i)}=\prod_{i=1}^k \frac{E_i(y_i')}{E_i(y_i)}D_i(x_i)$
$\bf{D}(\bf{x}|\bf{R})=\frac{\bf{E}(\bf{y'}|\bf{R}=\bf{I})}{\bf{E}(\bf{y}|\bf{R}=\bf{I})}\prod_{i=1}^k D_i(x_i)$
ずいぶんと簡単な式になったけれど、大丈夫だろうか…(未検算)
- 大丈夫そうなソース
- 正規分布・均一分布・指数分布・カイ二乗分布を例にそれらを組み合わせた複数の確率変数について実装してみた


MultiDistCorrelated<-function(k,xs,Dtypes,ArgsList,R){
	if(!is.matrix(xs))xs<-matrix(xs,nrow=1)
	# 特異値分解
	Re = sin(pi*R/2)
	Re.e=eigen(Re)
	if(min(Re.e[[1]])<0){
		print("Non-positive eigenvalue(s)!")
		print("Calculation was forced by non-zero eigenvalue(s) to 0.000000001.")
	}
	Re.e.ori=Re.e
	Re.e[[1]][Re.e[[1]]<0]=0.000000001
	# 標準正規分布の累積確率対応座標
	# Ds:確率密度
	ys.el<-xs
	Ds<-matrix(0,length(xs[,1]),k)
	for(i in 1:k){
		if(Dtypes[i]==1){
			#print(xs[,i])
			#print(pnorm(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))
			ys.el[,i]<-qnorm(pnorm(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))
			Ds[,i]<-dnorm(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2],log=TRUE)
		}else if(Dtypes[i]==2){
					#	print(xs[,i])

			#print(punif(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))
			ys.el[,i]<-qnorm(punif(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))
			Ds[,i]<-dunif(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2],log=TRUE)
		}else if(Dtypes[i]==3){
					#	print(xs[,i])

			#print(pexp(xs[,i],ArgsList[[i]][1]))
			ys.el[,i]<-qnorm(pexp(xs[,i],ArgsList[[i]][1]))
			Ds[,i]<-dexp(xs[,i],ArgsList[[i]][1],log=TRUE)
		}else if(Dtypes[i]==4){
				#		print(xs[,i])

			#print(pchisq(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))

			ys.el[,i]<-qnorm(pchisq(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2]))
			Ds[,i]<-dchisq(xs[,i],ArgsList[[i]][1],ArgsList[[i]][2],log=TRUE)
		}
	}
	#
	ys.sp<-t(Re.e[[2]]%*%diag(1/sqrt(Re.e[[1]]))%*%solve(Re.e[[2]])%*% t(ys.el))
	Pr.sp<-dnorm(ys.sp,log=TRUE)
	Pr.el<-dnorm(ys.el,log=TRUE)
	Pr<-exp(apply(Pr.sp-Pr.el+Ds,1,sum))
	return(list(Pr=Pr,Pr.sp=Pr.sp,Pr.el=Pr.el,Pr.xs=Ds))
}


# 任意の複数の確率分布を番号にてID化する
# 1 norm No.args=2 (args:mean var)
# 2 unif No.args=2 (args:min max)
# 3 exp  No.args=1 (args:lambda)
# 4 chisq No.args=2 (args:df ncp)
k<-2
Dtypes<-c(2,2)
ArgsList<-list(c(0,6),c(0,6),c(0,5),c(1,0))

x1<-seq(from=1,to=5,by=0.01)
x2<-seq(from=1,to=5,by=0.01)
x3<-1
x4<-1

xs<-as.matrix(expand.grid(x1,x2))
#x<-c(1,2,3,4)




R<-matrix(0.3,k,k)
diag(R)<-1
#R<-diag(rep(1,k))

outTestFx<-MultiDistCorrelated(k,xs,Dtypes,ArgsList,R)

filled.contour(matrix(outTestFx$Pr,length(x1),length(x2)))
min(outTestFx$Pr)-max(outTestFx$Pr)