R/シミュレーション

ランダムに正規分布を発生させる

 x = rnorm(10, mean=0, sd=1.5)

平均0、標準偏差1.5の正規分布のデータを10個発生させる。

乱数の種をセットする

set.seed(1)

指定した相関係数の2組のデータを発生させる

gendat2 <- function(nc, r)
{
# 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。
z <- matrix(rnorm(2*nc), ncol=2)
# 主成分分析を行い，主成分得点を求める。この時点で変数間の相関は完全に0である。
res <- eigen(r2 <- cor(z))
coeff <- solve(r2) %*% (sqrt(matrix(res$values, 2, 2, byrow=TRUE))*res$vectors)
z <- t((t(z)-apply(z, 2, mean))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff
# コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。
z %*% chol(matrix(c(1, r, r, 1), ncol=2))
}

使い方

Rho <- gendat2(10000, .5)

相関係数0.5の10000個のデータ対をランダムに作る。

相関係数行列で指定した乱数を発生させる

gendat <- function(     nc,                     # 標本サイズ
                        r)                      # 相関係数行列
{
        nv <- ncol(r)                           # 変数の個数
        z <- matrix(rnorm(nv*nc), ncol=nv)      # 仮のデータ行列を作る。この時点では変数間の相関は近似的に0である。
        r2 <- cor(z)                            # 相関係数行列
        res <- eigen(r2)                        # 主成分分析を行い，
        coeff <-  solve(r2) %*% (sqrt(matrix(res$values, nv, nv, byrow=TRUE))*res$vectors)
        z <- t((t(z)-colMeans(z))/sqrt(apply(z, 2, var)*(nc-1)/nc)) %*% coeff   #主成分得点を求める。この時点で変数間の相関は完全に0である。
        return(z %*% chol(r))                   # コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。
}

使用例

以下のような相関係数行列の３組のデータを２０組発生させる

10.5120.335
0.51210.467
0.3350.4671
> nc <- 20
> z <- gendat(nc, matrix(c(1, 0.512, 0.335, 0.512, 1, 0.467, 0.335, 0.467, 1), ncol=3))
> z

semパッケージの中に相関係数行列を指定する方法がある。

ガン <- read.moments(diag=TRUE,names=c("総熱量","肉類","乳製品","酒類",

"大腸ガン","直腸ガン"))

1.00
0.77 1.00
0.72 0.72 1.00
0.63 0.45 0.23 1.00
0.74 0.85 0.58 0.59 1.00
0.81 0.67 0.48 0.61 0.75 1.00

サンプリング

samplingに、10000の中から50個番号を選ぶ。 次に、Rhoというデータから、samplingで選んだ番号のデータを取り出す。

sampling <- sample(10000, 50)
temp <- Rho[sampling,]