# 多次元正規分布乱数を作る関数を作っておく library(MASS) library(lqmm) # 適当に対称行列を作って正定値行列を作る関数も作っておく # 多次元正規分布は中心座標と分散共分散行列とで決まる # 分散共分散行列のために正定値行列が必要 my.make.positive.definite <- function(n.v){# n.vは変数の個数 sigma <- matrix(rnorm(n.v^2),n.v,n.v) sigma <- sigma + t(sigma) sigma <- make.positive.definite(sigma) return(sigma) } # n.sampleは発生乱点数、n.vは変数の個数(次元) # mは中心座標、sigmaはn.v x n.v 正定値行列 my.mvrnorm <- function(n.sample,n.v,m=NULL,sigma=NULL){ if(is.null(m)){ m <- rnorm(n.v) } if(is.null(sigma)){ sigma <- my.make.positive.definite(n.v) } return(mvrnorm(n.sample,m,sigma)) } # サンプル数 n.s <- 1000 # 説明変数 ## 説明変数の数 n.v <- 5 ## 多次元正規分布とする M <- rnorm(n.v) # 重心 S <- my.make.positive.definite(n.v) # Sigma行列(分散共分散行列のようなもの) X <- my.mvrnorm(n.s,n.v,M,S) plot(as.data.frame(X)) library(rgl) tmp.X <- X[,1:3] tmp.X <- rbind(tmp.X,rep(max(tmp.X),3),rep(min(tmp.X),3)) # 3軸スケールをそろえるため plot3d(tmp.X) # sigma行列を適当に与えて、多次元正規分布(ただしゆがんでいる)を作った # 点の数をだんだん増やす n.ss <- c(10,25,50,75,100,250,500,750,1000,5000,10000) diff.M <- rep(0,length(n.ss)) for(i in 1:length(n.ss)){ n.s <- n.ss[i] X <- my.mvrnorm(n.s,n.v,M,S) # この乱点のラグビーボール状態を調べるのは、PCA # 標本重心 M.hat <- apply(X,2,mean) # 真値と比べる M.hat M pca.out <- prcomp(X) # 長軸からsdが並んで出る # sdを対角成分にした行列を作り K <- diag(pca.out$sdev) # 回転行列も出ているから、それで回す R <- pca.out$rotation P <- R %*% K Q <- P %*% t(P) # var-covとは各軸の内積のようなもの S diff.M[i] <- sum((Q-S)^2) } plot(log(n.ss),diff.M,type="l")
数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典