2013-12-03

相互に関連した乱数

n個の正規乱数が、平均はすべて0で、covariance matrix $\Sigma$ に従って得られるとき、地道に自分で作れば、適当な正の正方行列を作って、その固有値分解をして正規直交基底方向の拡大縮小成分と回転成分とに分けておき、多次元標準正規乱数から拡大縮小と回転とをさせることで作成できる
やってみよう

# 変数の数
n <- 4
library(GPArotation)
V <- Random.Start(n)
S <- diag(runif(n)*5)
sigma <- V%*%S%*%t(V)
#sigma <- matrix(c(4,2,2,3), ncol=2)
#sigma <- matrix(runif(n^2),ncol=n)
#sigma <- sigma+t(sigma)
# その固有値分解
eigen.out <- eigen(sigma)
eigen.out[[1]]
# 変数n個の正規乱数のベクトルをn.setセット作る
n.set <- 1000
X <- matrix(rnorm(n.set*n),ncol=n)
#plot(as.data.frame(X),xlim=range(X),ylim=range(X))
# 拡大縮小する
X. <- t(t(X) * sqrt(eigen.out[[1]]))
#plot(as.data.frame(X.),xlim=range(X.),ylim=range(X.))
# 回転する
X.. <- t(eigen.out[[2]] %*% t(X.))
plot(as.data.frame(X..),xlim=range(X..),ylim=range(X..))
# パッケージを読み込んでrmvnorm()を使ってみる
library(mvtnorm)
x <- rmvnorm(n=n.set, mean=rep(0,n), sigma=sigma)
colMeans(x)
colMeans(X..)
var(x)
var(X..)

# 別法
x <- rmvnorm(n=n.set, mean=rep(0,n), sigma=sigma, method="chol")
colMeans(x)
var(x)

このような乱数に、相互に独立した観測誤差を入れる
相互に共分散行列で関係付いた乱数は回転してやれば相互に独立した正規乱数として扱えるし、その回転において、観測誤差の間には関連が発生しないことを利用して、以下のようにして、観測結果から条件付期待値を算出してやることができる
やってみよう。上段は真値を横軸に、観測値・期待値とを縦軸にとった。下段は観測値を横軸に、真値・期待値を縦軸にとった。n個の乱数のうち3個を示したが、いずれも、真値~観測値の回帰直線と、条件付期待値との重なりがよいことが示されている

# 変数の数
n <- 4
library(GPArotation)
V <- Random.Start(n)
S <- diag(runif(n)*5)
sigma <- V%*%S%*%t(V)
#sigma <- matrix(c(4,2,2,3), ncol=2)
#sigma <- matrix(runif(n^2),ncol=n)
#sigma <- sigma+t(sigma)
# その固有値分解
eigen.out <- eigen(sigma)
eigen.out[[1]]
# 変数n個の正規乱数のベクトルをn.setセット作る
n.set <- 500
X <- matrix(rnorm(n.set*n),ncol=n)
#plot(as.data.frame(X),xlim=range(X),ylim=range(X))
# 拡大縮小する
X. <- t(t(X) * sqrt(eigen.out[[1]]))
#plot(as.data.frame(X.),xlim=range(X.),ylim=range(X.))
# 回転する
X.. <- t(eigen.out[[2]] %*% t(X.))
plot(as.data.frame(X..),xlim=range(X..),ylim=range(X..))

#########
# 上記の方法で発生させた共分散ありの多変数データをXに置きなおす
X <- X..
# 観測誤差は全変数で共通の正規乱数
s <- 1
Z <- rnorm(length(X..),0,s)
# 観測値
Y <- X + Z

# 観測値ありの下での条件付期待値は、共分散行列の固有値分解後の対角行列を
# 固有値そのままにするのは、「観測値そのものを期待値とする」ことに相当するが
# 条件付期待値の場合は、その分散を(真値の分散)/(真値の分散＋観測誤差の分散)ni
とすることなので、結局、固有値をその比で調整することが多変量化に相当する
X. <- t(eigen.out[[2]] %*% diag(eigen.out[[1]]/(eigen.out[[1]]+s^2)) %*% t(eigen.out[[2]]) %*% t(Y))

par(mfcol=c(2,3))
for(i in 1:3){


plot(X[,i],Y[,i],pch=20,xlab="真値",ylab="観測値・期待値")
points(X[,i],X.[,i],pch=20,col=2)
abline(0,1,col=3)
lm.out <- lm(Y[,i]~X[,i])
lm.out
abline(lm.out[[1]][1],lm.out[[1]][2],col=4)

plot(Y[,i],X[,i],pch=20,xlab="観測値",ylab="真値・期待値")
points(Y[,i],X.[,i],pch=20,col=2)
abline(0,1,col=3)

lm.out2 <- lm(X[,i]~Y[,i])
lm.out2
abline(lm.out2[[1]][1],lm.out2[[1]][2],col=4)
}

par(mfcol=c(1,1))