目的

特定の相関係数行列を持つ多変量データを生成する
（2変数データの場合は別のプログラムが便利）
MASSパッケージのmvrnorm関数でもできる

使用法

gendat(nc, r)

引数

nc	標本サイズ
r	相関係数行列

ソース

インストールは，以下の 1 行をコピーし，R コンソールにペーストする
source("http://aoki2.si.gunma-u.ac.jp/R/src/gendat.R", encoding="euc-jp")

# 特定の相関係数行列を持つ多変量データの生成
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))                   # コレスキー分解の結果をもとに，指定された相関係数行列 を持つように主成分得点を変換する。
}


使用例

> 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
       
        [,1]        [,2]       [,3]
   [1,]  2.21372763  0.90455236  1.5240454	# 生成されたデータ
   [2,] -0.45080094 -0.20986473  1.1174826
   [3,]  1.58487436 -1.06265744 -2.4285410
   [4,] -0.89568378 -1.13208582  0.3310756
   [5,]  1.36158887  0.84001818  1.7386909
   [6,]  0.39734615  1.04088405 -0.3364469
   [7,] -0.09416335 -0.20762585  1.0393830
   [8,] -1.22196036 -0.32578442 -0.6404094
   [9,] -0.84784564  0.08392657 -1.0298539
  [10,]  0.40092062  2.41337574  0.8425117
  [11,] -0.11916023  0.91027540 -0.6517597
  [12,] -1.27008358 -1.18811412 -0.2987666
  [13,] -2.18499030 -2.12059485 -1.6036881
  [14,]  0.64360886  0.01130282  0.1029740
  [15,]  0.34426143 -0.40514536  0.3413435
  [16,]  0.16459256  0.84058950  0.2445198
  [17,] -0.13113523  0.85887443 -0.2626515
  [18,]  0.33154003 -0.20170670 -0.4962152
  [19,] -0.30637310 -0.42684850  0.6376739
  [20,]  0.07973601 -0.62337127 -0.1713677

> apply(z, 2, mean)	# 平均値が 0 であることの確認
[1] -1.457168e-17 -8.326673e-17  1.387779e-18

> sqrt(apply(z, 2, var)*(nc-1)/nc)	# 標準偏差が 1 であることの確認
[1] 1 1 1

> cor(z)	# 指定した相関係数行列であることの確認
      [,1]  [,2]  [,3]
[1,] 1.000 0.512 0.335
[2,] 0.512 1.000 0.467
[3,] 0.335 0.467 1.000