相関のある２つの確率変数の同時分布 また続き

• こちらの続き
• 下の例は、指数分布の例
• これは、まだ、正規分布と指数分布との関係を使いきっていない状態
  • その点に対応するのは、明日の記事

# 指数分布の評価するべき値を指定
minX<-0
maxX<-5
kizami<-0.05

nV<-3 # 変数の数
# 変数ごとにパラメタ(期待値＝分散の平方根)を指定
af<-1:nV
af<-rep(1,nV)
# 評価値のリストを作る
tmpx<-seq(from=minX,to=maxX,by=kizami)
tmpx[1]<-0.001 # 0は定義域内だが、累積確率が0であって、対応する正規分布の値が取れないので少し大きくする
# 値リストについて変数の数の総順列作成
x<-list(tmpx)
X<-as.matrix(expand.grid(rep(x,nV)))
# 取った変数の値の累積確率を計算(指数分布)
pX<-X
for(i in 1:nV){
	pX[,i]<-pexp(X[,i],af[i])
}
# 累積確率が等しい標準正規分布の値
Xnorm<-qnorm(pX)

# 変数間の相関強度の行列
rmat<-matrix(0,nV,nV)
rmat[1,2]<-rmat[2,1]<-0.9
rmat[1,3]<-rmat[3,1]<-0.9
rmat[2,3]<-rmat[3,2]<-0.9
diag(rmat)<-1
# 固有値分解
ER = sin(pi*rmat/2)
ER.e=eigen(ER)
ER.e[[1]][ER.e[[1]]<0]=0
# 変数同士が独立だったとしたら、「非独立なXnorm」と同じ同時確率を持つ、標準正規分布の同時分布の値を算出
Xnorm2 = t((diag(1/sqrt(ER.e[[1]]))%*%t(solve(ER.e[[2]])))%*%t(Xnorm))
#Xnorm2 = Xnorm%*%diag(sqrt(ER.e[[1]]))%*%t(ER.e[[2]])
# 独立な同時分布の生起確率は掛け算なので計算可能
Pr<-exp(apply(dnorm(Xnorm2,log=TRUE),1,sum))
# 作業の過程で諸事情あって、XとTに入れ直し
T<-X
# 算出した生起確率を３変数次元のアレイに納める
PrArray<-array(Pr,c(length(tmpx),length(tmpx),length(tmpx)))
# 第１変数の値をアレイに納める
TArray1<-array(T[,1],c(length(tmpx),length(tmpx),length(tmpx)))
# 第２変数の値をアレイに納める
TArray2<-array(T[,2],c(length(tmpx),length(tmpx),length(tmpx)))

#persp(tmpx,tmpx,PrArray[,,1])
# プロットする
persp(TArray1[,1,4],TArray2[1,,4],PrArray[,,4],theta=80,phi=80,col="green")