【Rコード】多変量t分布のデータ発生【パッケージなし】

こんにちはたくまろです。今回はRコードで多変量t分布のデータ発生をパッケージなしで書いています。多変量正規分布のコードはこちらに載っています。多変量t分布を発生させる上で必要になります。

次元数が高いと計算速度が遅くなるので、変量t分布の定義に従って,Rコードを自分で書きました。仕組みとしては、多変量正規分布のデータを2つ発生させます。そのうち片方から、カイ2乗分布のデータを発生させます。データをカイ2乗分布のデータを自由度で割ったものを計算しています。

ここで、数学的に重要点ですが、t分布の共分散行列の定義が本によってことなります。ここでは、多変量t分布の共分散行列の値を入力値にしています。t分布を発生させる仮定で、正規分布を使うのですが、その正規分布の共分散行列を定義にしているものもあるので、どちらを基準にしているのか注意してください。

初めに、準備として多変量正規分布を発生させるこちらのRコードを読みこんでください。

datagene_norm <- function(mu,Sigma,N){
d <- dim(Sigma)[1]
Z <- matrix(0,d,N)
for (i in 1:N) {
Z[,i] <- rnorm(d,0,1)
}
gamma<- matrix(0,d,d)
diag(gamma) <- eigen(Sigma)$value^(1/2)
H <- eigen(Sigma)$vector
return(H%*%gamma%*%Z+mu)
}

こちらの関数を使って多変量t分布を発生させます。

###input###
# mu     mean vector
# Sigma  covariance matrix 
# N      sample size
# df(>2) degree of freedom
###output###
# N*d data matrix (N:sample, d:dimension)

datagene_t <- function(mu,Sigma,df,N){
if(df>2){
Sigma <- ((df-2)/df)*Sigma #correct covariance matrix
X <- datagene_norm(rep(0,d),Sigma,N) #generate N samples by N(0,Sigma)
Z <- datagene_norm(rep(0,df),diag(df),N) #generate N sample by N(0,I)
xi <- apply(Z^2,2,sum)/df #generate xi such that conformed Xi square distribution by Z
Y <- X%*%diag(sqrt(1/xi)) #make t by X and xi
t <- Y + mu #move in mu's direction
return(t)
}
else{print("Set df more than 3")}
}}

これを使って、多変量t分布のデータを発生させてみたいと思います。

d <- 100 #dimension 
N <- 10 #sample size 
df <- 10 #degree of freedom
mu <- rep(1,d) #mean vector 
Sigma <- matrix(0,d,d) #covariance matrix 
for(i in 1:d){ for(j in 1:d){ Sigma[i,j] <- 0.7^abs(i-j)  } }

X <- datagene_t(mu,Sigma,N) #Generate d*N data.
dim(X)
[1] 100 10

という風に出力され、正しい大きさになっています。

コメントを残す

メールアドレスが公開されることはありません。 * が付いている欄は必須項目です