2010-12-03

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

R 分布確率密度分布

こちらの続き
複数の確率変数がそれぞれ正規分布に従うときに、その同時分布は、変数間の相関の程度に応じて、「球」を「楕円」にすることで実現できる。それが、こちらで、正規分布だったらいいのにな、という話
球の楕円化は、行列でコレスキー分解とか固有値分解して正規直交基底の取り方を変えつつ、基底軸方向に伸び縮みさせること。いうれにしろ線形代数
以下の操作は、線形代数処理を使って、相互に相関のある正規分布変数に関して、多サンプルのデータを作った上で、その累積確率密度相当値を用いて、正規分布以外の確率変数に変換するというもの
例として指数分布でやってみたもの
10000サンプル、10変数、平均値が1:10 (したがって分散はその２乗の(1:10)^2)であって、変数相互の相関は相関係数として正規分布だったら0.7となるようにして発生させたもの
発生させた個々の変数の平均と分散は、指定したものになっていることを確認し
個々の変数の分布が指数分布様であることを１変数の度数分布で確認し
２変数の同時分布を２次元ヒストグラムとして描いてある
結果として得られた相関係数はだいたい、0.87くらい
さらに、変数間の相関の程度をばらばらにしてやっても、まずまず、うまく相関がとれる

mvrexp = function(n=10000, m=20, r=0.1, method="eigen", af=NULL){
rho<-NULL
ER<-NULL
 # expected correlation matrix
  if(is.matrix(r)){
    ER = sin(pi*r/2)
    #m=dim(ER)[1]
  }else{
    ER = matrix(rep(sin(pi*r/2), m^2), m,m)
    diag(ER) = 1
  }
  # n m-variate standerd normal random numbers
  X = matrix(rnorm(n*m),n)
  # m-variate normal random numbers with given correlation
  if(method=="cholesky"){
		X = X%*%chol(ER)
  }else if(method=="eigen"){
		ER.e=eigen(ER)
		ER.e[[1]][ER.e[[1]]<0]=0
		X = X%*%diag(sqrt(ER.e[[1]]))%*%t(ER.e[[2]])
  }
 
  # transformation to other variables
  if(is.null(af)){
		af = rep(1, m)
  }else if(length(af)==1){
		af = rep(af, m)
  }
  for(i in 1:m){
	#X[,i] = as.numeric(X[,i]<=qnorm(af[i]))
	X[,i] = qexp(pnorm(X[,i]),af[i])
	#X[,i]<-qunif(pnorm(X[,i]))
  }
return(X)
}
n<-10000
m<-10
af<-1:m
r<-0.7
# correlation matrixは実データから作らないと、「そんなのありえない！」になるので
tmpdata<-matrix(runif(n*m),n,m)
for(i in 1:m){
 tmpdata[,i]<-tmpdata[,i]+tmpdata[,sample(1:m,1)]*rpois(1,4)
}
rmat<-cor(tmpdata)
X<-mvrexp(n=n,m=m,af=1/af,r=rmat)
#par(mfcol=c(2,4))
layout(matrix(c(1,2,3,4,5,6,6,7,7,7,6,6,7,7,7),3,5,byrow=TRUE))
#image(X)
cor(X)
image(rmat)
image(cor(X))
plot(c(cor(X)),c(rmat))
# mean of m variables
apply(X,2,mean)
# var of m variables
apply(X,2,var)
#plot(as.data.frame(cbind(af,apply(X,2,mean),sqrt(apply(X,2,var)))))
plot(af,apply(X,2,mean))
plot(apply(X,2,mean),sqrt(apply(X,2,var)))
hist(X[,1])
library(gregmisc)
h2d <- hist2d(X[,1],X[,2], show=FALSE,same.scale=TRUE,nbins=c(10,10))
#filled.contour( h2d$x, h2d$y, h2d$counts, nlevels=9,col=gray((8:0)/8) ) # 2 次元ヒストグラムを濃淡で
persp( h2d$x, h2d$y, h2d$counts,ticktype="detailed",theta=50, phi=30,shade=0.5, col="red")
#plot(X[,1],X[,2])
par(mfcol=c(1,1))
#par(def.par)
#plot(c(cor(X)),c(rmat))