ぱらぱらめくる『Rによる計算機統計学』と『生物学のための計算機統計学』
- 作者: Maria L..Rizzo,石井一夫,村田真樹
- 出版社/メーカー: オーム社
- 発売日: 2011/05/26
- メディア: 単行本（ソフトカバー）
- 購入: 1人 クリック: 36回
- この商品を含むブログ (5件) を見る
- 『Rによる計算機統計学』の著者のサイト
- そこから辿るグラフィクス一覧サイト
- 本の内容など
- その本の中のソースのダウンロードサイト
- 第３章 確率変数の発生
- 逆変換法(関数と逆関数を使う)
- 採択・棄却法
- 変換法(確率分布間の関係を使う(こちら(とその元目次(こちら))からリンクを張ったこちらのクリッカブル・チャートを))
- 総和・畳み込みと混合
- 分布の混合に関するikuro's homepage 記事はこちら
n<-1000 m1<-0 m2<-10 sd1<-1 sd2<-5 x1<-rnorm(n,m1,sd1) x2<-rnorm(n,m2,sd2) tatamikomi<-x1+x2 k<-as.integer(runif(n)>0.5) kongou<-x1*k+x2*(1-k) par(mfcol=c(1,2)) hist(tatamikomi,prob=TRUE) hist(kongou,prob=TRUE) par(mfcol=c(1,1))
-
- 多変量分布
- 多変量正規乱数
- 多変量分布
### Example 3.18 (Choleski factorization method) # 標本データから各確率変数の平均と分散共分散行列を出して、それをもとに多変量正規分布乱数を発生する rmvn.Choleski <- function(n, mu, Sigma) { # generate n random vectors from MVN(mu, Sigma) # dimension is inferred from mu and Sigma d <- length(mu) Q <- chol(Sigma) # Choleski factorization of Sigma Z <- matrix(rnorm(n*d), nrow=n, ncol=d) X <- Z %*% Q + matrix(mu, n, d, byrow=TRUE) X } #generating the samples according to the mean and covariance #structure as the four-dimensional iris virginica data y <- subset(x=iris, Species=="virginica")[, 1:4] mu <- colMeans(y) Sigma <- cov(y) mu Sigma #now generate MVN data with this mean and covariance X <- rmvn.Choleski(200, mu, Sigma) pairs(X)
-
-
-
- こちらでMASSパッケージのmvrnorm()を紹介されているので、その中身と比較してみよう
- eigen()を使うかchol()を使うかは、ひとまず置いて置くとして、
- 教科書の方は『どういうつもりでこのソースを書いていますよ』という書き方になっていて、
- mvrnorm()の方はRのパッケージにするにあたっての不備がないようにする行が多い
- こちらでMASSパッケージのmvrnorm()を紹介されているので、その中身と比較してみよう
-
-
library(MASS) mvrnorm
> mvrnorm
function (n = 1, mu, Sigma, tol = 1e-06, empirical = FALSE)
{
p <- length(mu)
if (!all(dim(Sigma) == c(p, p)))
stop("incompatible arguments")
eS <- eigen(Sigma, symmetric = TRUE, EISPACK = TRUE)
ev <- eS$values
if (!all(ev >= -tol * abs(ev[1L])))
stop("'Sigma' is not positive definite")
X <- matrix(rnorm(p * n), n)
if (empirical) {
X <- scale(X, TRUE, FALSE)
X <- X %*% svd(X, nu = 0)$v
X <- scale(X, FALSE, TRUE)
}
X <- drop(mu) + eS$vectors %*% diag(sqrt(pmax(ev, 0)), p) %*%
t(X)
nm <- names(mu)
if (is.null(nm) && !is.null(dn <- dimnames(Sigma)))
nm <- dn[[1L]]
dimnames(X) <- list(nm, NULL)
if (n == 1)
drop(X)
else t(X)
}
<bytecode: 0x09bde874>
<environment: namespace:MASS>
-
-
-
- と思ったら、本では、（rmvn.eigen(),rmvn.svd()も作って、それらとともに）パフォーマンスを比較していた
-
-
> set.seed(100) > system.time(for (i in 1:N) + rmvn.eigen(n, mu, cov(matrix(rnorm(n*d), n, d)))) ユーザ システム 経過 2.95 0.00 2.95 > set.seed(100) > system.time(for (i in 1:N) + rmvn.svd(n, mu, cov(matrix(rnorm(n*d), n, d)))) ユーザ システム 経過 3.79 0.00 3.81 > set.seed(100) > system.time(for (i in 1:N) + rmvn.Choleski(n, mu, cov(matrix(rnorm(n*d), n, d)))) ユーザ システム 経過 2.2 0.0 2.2 > set.seed(100) > system.time(for (i in 1:N) + mvrnorm(n, mu, cov(matrix(rnorm(n*d), n, d)))) ユーザ システム 経過 3.03 0.00 3.02 > set.seed(100) > system.time(for (i in 1:N) + rmvnorm(n, mu, cov(matrix(rnorm(n*d), n, d)))) ユーザ システム 経過 3.82 0.00 3.84 > set.seed(100) > system.time(for (i in 1:N) + cov(matrix(rnorm(n*d), n, d))) ユーザ システム 経過 0.97 0.00 0.97
m<-10000 x<-runif(m) theta.hat<-mean(exp(-x)) theta.hat 1-exp(-1)
> theta.hat [1] 0.6321605 > 1-exp(-1) [1] 0.6321206
m<-10000 L<-2 U<-4 x<-runif(m,min=L,max=U) theta.hat<-mean(exp(-x))*(U-L)
> theta.hat [1] 0.1164698 > exp(-L)-exp(-U) [1] 0.1170196
- 章横断的に「モンテカルロ法について」のコメント
- 所定の確率モデルから、パラメトリックなブートストラップとも呼ばれる反復サンプリング
- 観測標本からのリサンプリングに基づくノンパラメトリックなブートストラップ法
- マルコフ連鎖モンテカルロ法
- 第６章 推論におけるモンテカルロ法
- 第７章 ブートストラップ法・ジャックナイフ法
- リサンプリング
- ノンパラメトリックなブートストラップ法
- 第８章 並べ替え検定
- リサンプリング２
- 第９章 マルコフ連鎖モンテカルロ法
- 第１０章 確率密度推定
- 第１１章 Ｒによる数値解析はこちら
生物学のための計算統計学 ?最尤法，ブートストラップ，無作為化法?
- 作者: Derek A. Roff,野間口眞太郎
- 出版社/メーカー: 共立出版
- 発売日: 2011/03/10
- メディア: 単行本
- クリック: 24回
- この商品を含むブログ (8件) を見る
