乱数をつかって任意のデータセットを作るなど

Rによるやさしい統計学を参考に，乱数を使って任意のデータセットを作る関数を作ってみた。

まずは基本。任意の相関係数を持つ二変数を生成する関数。

#任意の相関係数を持つ二変数を生成する関数
cor2val <- function(rho,size){
  X <- rnorm(size)
  e1    <- rnorm(size)
  e2    <- rnorm(size)
  V1 <- sqrt(rho)*X+sqrt(1-rho)*e1
  V2 <- sqrt(rho)*X+sqrt(1-rho)*e2
  
  result <- data.frame(cbind(V1,V2))
  return(result)
}

使い方は簡単で，相関係数と作るデータセットのサイズを引数で渡すだけ。

obs <- cor2val(size=500,rho=0.3)
summary(obs)
cor(obs)

次は任意の相関行列をもつデータセットを作る関数。
#任意の相関行列をもつデータセットを作る関数

corMatVal <- function(size,n.vals,cor.mat,corvec=FALSE){
  Rmat <- matrix(rnorm(n=size*n.vals),nrow=size)
  #if given cor.mat is just vector,make correlational matrix
  if(corvec==TRUE){
    temp.mat <- diag(n.vals)
    temp.mat[(upper.tri)(temp.mat,diag=FALSE)] <- cor.mat
    temp.cor <- (t(temp.mat)+temp.mat)
    diag(temp.cor) <- 1
  }else{
    temp.cor <- cor.mat
  }
  Upper <- chol(temp.cor)
  result <- as.data.frame(Rmat %*% Upper)
  return(result)
}

引数に生成するデータサイズ，変数の数，目標となる相関行列を渡す。

target.Mat <- matrix(c(1.0,0.4,0.3,0.2,0.3,
                       0.4,1.0,0.8,0.1,0.2,
                       0.3,0.8,1.0,0.4,0.5,
                       0.2,0.1,0.4,1.0,0.1,
                       0.3,0.2,0.5,0.1,1.0),nrow=n.vals)
obs2 <- corMatVal(500,5,target.Mat)

ここでちょっと一工夫した。相関行列は対称で対角が1に決まってるんだから，全部入れるのめんどくさい。例えば3変数の相関行列なら，r_12,r_13,r_23の情報だけでいいじゃない。ということで，そういう横着をしたい人は，目標行列のところに額各要素のベクトルを渡して，corvecオプションをTRUEにすればいいようにした。

obs3 <- corMatVal(500,3,c(0.2,0.3,0.4),corvec=TRUE)

最後に，任意の因子負荷行列を持つデータセットの生成関数。

#任意の因子負荷行列をもつデータセットを作る関数

FloadVal <- function(size,load.mat){
  n.factors <- ncol(load.mat)
  n.vals <- nrow(load.mat)
  uni.vec <- apply(load.mat,1,function(x) 1-sum(x^2))
  uni.mat <- diag(sqrt(uni.vec))
  fscore <- matrix(rnorm(n.factors*size),nrow=n.factors)
  uniscore <- matrix(rnorm(n.vals*size),nrow=n.vals)
  obs <- t(load.mat%*%fscore+uni.mat%*%uniscore)
  result <- as.data.frame(obs)
  return(result)
}

こんな風にして使う。

load.mat <- matrix(c(0.6,0.6,0.2,0.2,0.3,0.3,0.4,0.4,0.5,0.5),nrow=5)
obs4 <- FloadVal(size=50000,load.mat)

備え付けの因子分析関数で確認。バッチリ。

factanal(obs4,2)

もっとも，因子間相関がない仮定。実際に斜交回転で分析すると相関が出て負荷量もおかしくなっちゃう。
これについてはちょっと今度手を入れないと行けないですなー。