「統計学関連なんでもあり」の過去ログ--- 042

No.08412　データの作成　　【空星】　2008/11/27(Thu) 13:40

指定した相関係数行列，平均，標準偏差の下で項目数＞データ数のデータを作成したいのですが，例えばデータ数10，項目数20のデータを作成する場合，

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

convert <- function(z, m, s) return(t(t(z)*s+m)) # 標準化の逆

nc<-10
z<-gendat(nc,matrix(c(
1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1
),ncol=20))
m <- c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0) # 必要とする平均値
s <- c(1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1,1)  # 必要とする標準偏差
z2 <- convert(z, m, s) # 標準化の逆を行う
colMeans(z2) # 平均値の確認
sd(z2)*sqrt((nc-1)/nc) # 標準偏差の確認

でデータを作成しようと思っているのですが
「以下にエラーsolve.default(r2) : システムは数値的に特異です：条件数の逆数 = 1.16668e-018」
とエラーがでてしまいうまくいきません。どうすれば，項目数＞データ数のデータを作れるのでしょうか？ご教授いただければ助かります。
宜しくお願いします。

No.08413　Re: データの作成　　【青木繁伸】　2008/11/27(Thu) 14:32

相関係数が単位行列（つまり，どの変数も他の変数との相関は0）となるようなデータを求めたいのですか？

そうであるとして，変数が20個で，所定の相関係数行列（単位行列でなくても）を持つデータを定めるためには標本サイズは21以上でなくてはなりません。「システムは数値的に特異です」というRのメッセージは不適切というか言葉遣いに難があるのですが，要するに計算の途中で特異行列が出てきてしまって計算を継続することができないということです。特異行列とは要するに逆行列が計算できない（一次従属な行・列を含む）ということです）。ですから上のようなことを行うためには
> convert <- function(z, m, s) return(t(t(z)*s+m)) # 標準化の逆
> nc <- 21
> z <- gendat(nc, diag(20)) # 単位行列は diag 関数で作る（間違いが無く，簡単）
> colMeans(z)
 [1] -2.114711e-17  2.643388e-18  8.591012e-18  2.114711e-17  3.886607e-17
 [6] -2.841642e-17  3.634659e-17 -9.086647e-18  4.493760e-17 -5.286776e-18
[11]  4.229421e-17  6.608470e-18  1.321694e-17  5.286776e-17 -2.907727e-17
[16] -1.850372e-17 -2.521131e-16 -9.119689e-17 -2.590520e-16  1.321694e-17
> sd(z)*sqrt((nc-1)/nc)
 [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
> m <- 1:20*0.1
> s <- 1:20*0.05
> z2 <- convert(z, m, s)
> colMeans(z2)
 [1] 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1.0 1.1 1.2 1.3 1.4 1.5 1.6 1.7 1.8 1.9
[20] 2.0
> sd(z2)*sqrt((nc-1)/nc)
 [1] 0.05 0.10 0.15 0.20 0.25 0.30 0.35 0.40 0.45 0.50 0.55 0.60 0.65 0.70 0.75
[16] 0.80 0.85 0.90 0.95 1.00
のようになります。なお，これと同じことはMASSパッケージのmvrnorm関数でもできます（というか，mvrnormを使う方が若干簡単）。