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