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

No.04393　データの生成　　【空星】　2007/09/24(Mon) 11:39

以前，指定した相関係数行列のもとで，指定した平均値と標準偏差のデータの生成をExcelのgendatのファイルを教えてもらい解決したのですが，このデータの発生をRを使って行いたいのですが，青木先生のRによる統計処理のデータ生成・変数変換を参考にしたのですが，うまくいきません。このようなソースはあるのでしょうか？よろしくお願いします！

No.04396　Re: データの生成　　【青木繁伸】　2007/09/24(Mon) 18:46

どのようにうまくいかないのですか?

＞このようなソースはあるのでしょうか

あそこにあるやつは，全部ソース付きですよ

No.04399　Re: データの生成　　【空星】　2007/09/25(Tue) 09:52

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,0,0,1,0,0,0,1),ncol=3))

このソースで，matrix(c(1,0,0,0,1,0,0,0,1)の相関係数行列で，平均が0，標準偏差が1のデータ数が20で3列のデータを生成してから，3列目を平均3で標準偏差1にする場合，

gendat2<-function(mu=0,sigma=1){
return((z[,3]-mean(z[,3]))/sd(z[,3])*sigma+mu)
}

z[,3]<-gendat2(3,1)

とするのですが，gendat2の関数の中身をz[,3]のように書き換えているのですが，指定した平均と標準偏差にする度に書き換えないでできないのかなと思っていまして。列が1つの場合を変えるのはいいのですが，これだと2つの列の平均，標準偏差を変えると相関係数行列に影響してくるのではないかと思っていまして。アドバイスよろしくお願いします！

No.04400　Re: データの生成　　【青木繁伸】　2007/09/25(Tue) 10:41
> gendat2<-function(mu=0,sigma=1){
> return((z[,3]-mean(z[,3]))/sd(z[,3])*sigma+mu)
> }
z は完全に標準化されているので (z[,3]-mean(z[,3]))/sd(z[,3]) は不要です。return(z[,3]*sigma+mu) でよいです。
なお，z の分散は不偏分散ではない方ですので，sd(z) をしても1にはなりません。

> 指定した平均と標準偏差にする度に書き換えないでできないのかな

希望する平均値ベクトル，標準偏差ベクトルを m, s としたとき t(t(z)*s+m) で望むものはできるでしょう。

毎回このようなものを書くのがいやだと言うことならば，

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

みたいなものを書く（gendat() に m, s を渡してやっても良いのだが）

> 列が1つの場合を変えるのはいいのですが，これだと2つの列の平均，標準偏差を変えると相関係数行列に影響してくるのではないか

各変数をどのように線形変換しようと，相関係数行列は不変です
> nc<-20
> z<-gendat(nc,matrix(c(1,0,0,0,1,0,0,0,1),ncol=3))
> m <- c(50, 60, 90) # 必要とする平均値
> s <- c(10, 20, 30)  # 必要とする標準偏差
> z2 <- convert(z, m, s) # 標準化の逆を行う
> colMeans(z2) # 平均値の確認
[1] 50 60 90
> sd(z2)*sqrt((nc-1)/nc) # 標準偏差の確認
[1] 10 20 30
> cor(z2) # 相関係数の確認
              [,1]          [,2]         [,3]
[1,]  1.000000e+00 -2.713940e-16 4.409436e-16
[2,] -2.713940e-16  1.000000e+00 5.344752e-16
[3,]  4.409436e-16  5.344752e-16 1.000000e+00

No.04422　Re: データの生成　　【空星】　2007/09/27(Thu) 15:44

無事うまくいくことができました。どうもありがとうございました。