# 同一染色体上にある多型の数Ns Ns<-4 # Ns個のマーカーのそれぞれのアレルの数をNaS NaS<-c(3,4,5,2) # アレル頻度をPaS PaS<-list() for(i in 1:Ns){ PaSi<-c(rdirichlet(1,rep(1,NaS[i]))) } NaS PaS hap1<-hap2<-rep(0,Ns) for(i in 1:Ns){ hap1[i]<-sa…

> legend(x=par()$usr[2],y=par()$usr[4],legend=colnames(dat),pch=22,lty=0,xjust=-.1+ ,pt.bg=cm.colors(3),col=1:3,pt.cex=2,x.intersp=2,y.intersp=1.5,title="xpd=TRUE")プログラムの2行目以降（継続行）が "," で始まるというスタイルは，見たことが…

figure チャンクを使って図を挿入するのは記述が面倒くさい。そこで，以下のような関数を作っておき，results=tex を使って作図した結果を挿入する。 xfigure <- function(src, fn, caption="", label=fn, size=320, width=6, height=6, mgp=c(1.8, 0.6, 0),…

if 文 というタイトルのついたページ 例を挙げているのだから，揚げ足取ってもしようがないのだけど，パズルとして 以下のようなものは， > set.seed(123) > x <- sample(9,12,replace=TRUE) > x <- matrix(x,nrow=2,dimnames=list(1:2,letters[1:6])) > for…

変わった散布図 symbol 関数，そして引数 pch について少し にて > x <- cbind("a"=1:5,"b"=5:1,"c"=6:10,"d"=10:6,"e"=sample(10,5)/10) > symbols(x[,1:2],circles=x[,1]/10,inches=FALSE,main="circles")散布図において，半径（直径）が度数を表すように…

棒グラフのいろいろ・その2 にて ここに示されているようなグラフを，ダイナマイトプロットと呼ぶ。描いてはいけないグラフの一つ。 なお，グラフを jpeg デバイスで描くのもどうかと思う。せめて png で描こう。 また，プログラムで，カンマの後に空白がな…

ヒストグラムと箱髭図 にて > truehist関数は総面積が1に基準化されるので、縦軸が異なります。 というのはうそ。freq=FALSE にすれば同じこと。 一番違うのは，hist では right=TRUE がデフォルトになっていて，一般常識と異なる点（事実，リンク先のページ…

image関数でイメージする(R Advent Calenderの番がきた) にて > (x <- matrix(1:12,ncol=3)) [,1] [,2] [,3][1,] 1 5 9[2,] 2 6 10[3,] 3 7 11[4,] 4 8 12> (y <- t(x[nrow(x):1,ncol(x):1])[ncol(x):1,]) [,1] [,2] [,3] [,4][1,] 4 3 2 1[2,] 8 7 6 5[3,] …

Sweave-12 で書いた ae フォントを使うかどうかの件 やはり，困るという人もいるようで，マニュアルにちゃんと記載がある。 ‘Sweave.sty’ also supports the [noae] option, which suppresses the use of the ae package, the use of which may interfere wi…

「マクネマー検定についてのページ」で epibasix パッケージに mcNemar 関数があるのを知った。しかし，この関数の実装は余り良くない。 (1) P 値が 0.001 より小さい場合， McNemar's Chi^2 Statistic (corrected for continuity) = 12.345 which has a p-v…

Sweave で作った pdf と，LaTeX だけで作った pdf を比べて，前者と後者のフォントの大きさ（当然ながら，版面の大きさとか行当たりの文字数など）が異なることに気づいた。LaTeX のソースに \usepackage{Sweave} を加えるだけで，変になることがわかった。 …

apply関数は遅かった・・・ 確かに apply は遅い。その代わりに colMeans, rowSums を使って定義通りに書くと，ヘタに行列演算するより速い。colMeans, rowSums を使って書いたものと，結果の格納領域を事前に確保して行列演算をするのは，ほとんど同じくら…

今回は，R プログラミングではなく，統計学の方でダメ出し例題をやってみる その2 で，> ここでは、抗生物質を使うとプラセボより感染率は小さくなるはずだ、という期待が込められているので片側検定を行なっていると思われる。のに> （ここではカイ二乗検定…

相変わらず，「空白などという無駄なものは置くものか」という，読む気のしないグッチャリプログラム 元のプログラムは 98 秒かかる（作者の環境では 3 分ぐらいかかるそうだ）が， pvalue <- sum(apply(boot, 1, mean) > 37)/l を pvalue <- sum(rowMeans(b…

int64 パッケージ > library(int64)> a <- as.int64(1)> for (i in 1:300) {a <- a*i; print(a)}[1] 1[1] 2[1] 6[1] 24[1] 120[1] 720[1] 5040[1] 40320[1] 362880[1] 3628800[1] 39916800[1] 479001600[1] 6227020800[1] 87178291200[1] 1307674368000[1] 2…

等式を成り立たせようさっきのと同じように解ける。ユニークな解は1つだけど。 > library(e1071)> invisible(apply(permutations(9), 1, function(i) {+ st <- paste(c("", "/", "", "+", "/", "", "+", "/", ""), i, collapse="", sep="")+ if (eval(parse(…

式を成り立たせる問題。簡単すぎる。×を*，÷を/，=を==とする。 > library(e1071)> op <- c("+", "-", "*", "/", "==")> invisible(apply(permutations(5), 1, function(i) {+ st <- paste(paste(5:1, op[i], collapse="", sep=""), "0", sep="")+ if (eval(…

R 2.14.0 から，pbirthday() and qbirthday() now use exact calculations for coincident = 2.になった。 R 2.13.1 では，pbirthday() and qbirthday() did not implement the algorithm exactly as given in their reference and so were unnecessarily in…

既約分数クイズ にあるもの。 0/q, q/q は自明なので，除く。いかにもRらしい解??。ただし，メモリ節約も最適化も何にも考えない。既存関数を利用し，ひたすら短く難解なプログラムを目指す。 > func <- function(d0)+ {+ euclid <- function(m, n)+ {+ m0 <…

なぜ廃止されたのかな？ストイックな理由なんだろうか。 現在の sd は > sdfunction (x, na.rm = FALSE) { if (is.matrix(x)) { msg <- "sd(<matrix>) is deprecated.\n Use apply(*, 2, sd) instead." warning(paste(msg, collapse = ""), call. = FALSE, domain = </matrix>…

sd 関数は，もはや，行列やデータフレームを対象としない。 > class(d) [1] "matrix" > sd(d) [1] 1 1 警告メッセージ： sd() is deprecated. Use apply(*, 2, sd) instead. > sd(data.frame(d)) X1 X2 1 1 警告メッセージ： sd() is deprecated. Use sapply…

> x <- data.frame(x=1:5, y=c(3,2,5,4,3))> sapply(x, factor) x y [1,] "1" "3"[2,] "2" "2"[3,] "3" "5"[4,] "4" "4"[5,] "5" "3"

平均値を求めるくらいなら，何の差もないだろうと... 最初に，素直にベクトル要素に添え字でアクセスしてみたら，mean 関数の速度の半分ほど（要するに，遅いと言うこと）。 「これではならじ」と，ポインタを使って書いたら， src <- ' Rcpp::NumericVector…

e1071 に permutations という関数がある。実装はいろいろあるけど，簡単そうな奥村先生に典拠するアルゴリズムを使って書いてみた。ただし，n! なんて，すぐにメモリ制限を超えてしまうので，実際の利用はたかだか n <= 10 しかないので，どんなに速いアル…

2011-11-06 ２次計画問題 http://d.hatena.ne.jp/ryamada22/20111106/1320544702 segments は，ベクトル化されている# for(i in 1:N){# segments(i,D[i,1],i,D[i,2],col=3)#}ではなくsegments(1:N, D[, 1], 1:N, D[, 2], col=3)しかし，for(i in 1:N){ for(j…

R 2.13のコンパイラをつかってみる(3)http://ito-hi.blog.so-net.ne.jp/2011-07-09-2> set.seed(1234)> n <- 10000000 # 速すぎて速度が測れないので多めに設定> > x <- rnorm(n, 0, 1)> > # ctest2 は test2 をコンパイルしたもの> test2 <- function(x) {+…

R 2.13.0のコンパイラをつかってみる (2) Cλhttp://ito-hi.blog.so-net.ne.jp/2011-04-21Morisita の Cλ を求めるのに，行ごとの組合せで n.cmm × n.cmm 回 clambda 関数を呼んでいるが，無駄。対称行列なのだから，せめて半分だけ計算する。以下，原著者に…

移動平均は stats::filter で行える。 これを，inline で書いて見ようというのが今回のお題。 NA は NA_REAL を使う（IntegerVector なら，NA_INTEGER。他に，Inf, -Inf は R_PosInf, R_NegInf。） ポインターは double *pz; ... pz = REAL(z); のように使う…

コメントするのにログインまたは登録しろと。。。そんなのやだね。登録者数を増やしたいのか。 http://d.hatena.ne.jp/tsutatsutatsuta/20111024 行列の各要素ごとに異なる引数を与えて，同じ関数を適用 applyを使って，行列の行や列ごとに，異なる引数を与…

簡単に使えて，効果もそれなりに見込める > func <- function(xa, xb) {+ xab <- integer(length(c(xa, xb))-1)+ for (i in seq.int(xa))+ for (j in seq.int(xb))+ xab[i+j-1] <- xab[i+j-1]+xa[i]*xb[j]+ return(xab)+ }> system.time(a <- func(1:100, 10…