関数を作ってみました。 HS.tapply <- function(x, index, type) { # create index indicate matrix index <- as.numeric(factor(index)) mindex <- NULL for(i in 1:length(unique(index))) { mindex <- cbind(mindex, index-(i-1)) } if (type == 5 || typ…

うーん、こりゃ参ったねｗ n <- 20 rep <- 100000 d <- matrix(rnorm(n*rep), nrow=n, ncol=rep ) # apply Rprof() varA <- apply(d, 2, var) Rprof(NULL) print(a <- summaryRprof()$sampling.time) # matrix operation (1) Rprof() varB <- NULL for(i in …

前回書いた結果をもとにして、ちゃんと近似できているのか試してみることにしました。 今回は Welch (1938) の (2) 式の統計量 \[ v=\frac{\bar{X}_1-\bar{X}_2}{\sqrt{\frac{V_1}{n_1}+\frac{V_2}{n_2}}} \] を導出する過程で出てくる分布の近似を考えてみ…

シミュレーション中に進捗状況を監視。 path <- "C:/"; set.seed(091021) r <- 20 m <- c(0, 0, 0.3) s <- c(.1, .1, .1) n <- c(10, 10, 10) d <- do.call("rbind", mapply(function(n, m, s) matrix(rnorm(r * n, mean = m, sd = s), nrow=n, ncol=r), n, …

Rでシミュレーションを行おうと思って、「群変数が共通で、データを乱数として r セット生成する」スクリプトをループを使わずに書こうと思ったら、なかなかうまくいかない。 ↑を逆(「データが共通で、群変数を乱数として r セット生成する」)にするとリサン…

通常フォントファミリーで指定してしまうので、フォントを混在させることはあまりありません。 mtext関数を駆使する事で、複数のフォントを使ってグラフを作成することができます。 postscript出力では未チェックですが、Windows環境でemf出力するか、ビット…

expression()関数中で書体の切り替えができたらしい、何で今まで気づかなかったんだろう・・・。 plain(x) draw x in normal font bold(x) draw x in bold font italic(x) draw x in italic font bolditalic(x) draw x in bolditalic font 他の基本的な数学…

正しいのかどうか分かりませんが、一応できたのでメモっておきます。 URW Garamondを使ってみよう。 まず、Rをインストールしたディレクトリの中の /R/R-X.X.X/library/grDevices/afm に使いたいフォントのafmファイルを一つずつgzファイルにしてからコピー…

barplotでやると盛大にはみ出ます・・・ Ety <- function(y) { n <- 20 return( (y+1)/(n+2) ) } par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c(2.0, 0.7, 0)) par(bty="l") names <- replace(seq(0, 20), which(seq(0, 20) %% 5 != 0), NA) barplot( Ety(seq(0,20)…

plot()→abline()などで参照線を引くと、参照線のはずなのにグラフの線より前に表示されてしまう。 なので、panel.firstオプションを使うと、参照線→グラフ線の順に描画してくれる。 xx <- runif(100) yy <- runif(100) par(mar=c(3.2, 3.2, 1.5, 0.5), mgp=c…

Rのデフォルトだとちょっと軸ラベルが離れすぎているし、余白も多めなので調べてみた。 グラフィックスパラメータの mpg で、軸ラベル・軸目盛・軸線のバランスを整える事ができる様だ。 ヘルプには以下のような感じに記述してあった。 軸からラベル、軸から…

Rでもバッチ実行ディレクトリにファイルを保存したかったので調べてみた。 Rscriptを使うと出来るようだ。 Rscript [--options] [-e expr] inputfile [args]のように、argsを引数として渡すことができる。 Rscript -e "commandArgs(); proc.time()" arg1 arg…

Windowsにインストールされているフォントを適当に指定してfamilyとして定義すると使えるようになるみたい。 windowsFonts(Times="Times New Roman") par(family="Times") 定義済みのフォントは windowsFonts() で表示される。 n <- 30 y <- 2 Bin <- functi…

不完全ベータ関数の近似はまとめるほど理解できなかったので、また時間があるときにやろう。 今日はラプラス近似が何をしてるかを少しグラフィカルに検討した。 \[ \int_{\mathbf b\in R^q}\exp\{h(\mathbf b)\}d\mathbf b \simeq \exp\{h(\hat{\mathbf b})\…

メモが目的なので要点だけ。 dev.copy2eps()だと何故かうまくいかなかったので、postscript()を使った。 postscript( file="c:/test.eps", width = 9, height = 9, horizontal = FALSE, onefile = FALSE, paper = "special", family="AvantGarde" ) family="…

2〜20次元関数のAdaptive Quadrature(適応的求積法)ができるパッケージ。 library(adapt) # Example of p - dimensional normal distribution mu <- rbind(c(0, 0, 0)) sigma <- rbind(c(1, 0, 0), c(0, 1, 0), c(0, 0, 1)) pdfnm <- function(z) { 1 / ( (2…