ダメ出し：ryamada本のp.213の図12.3をRで描く

どうも，この一派（？）は，プログラムはスペースを使わずギッチリ詰めて書くとか，par(new=T) を使って（せめて par(new=TRUE) とするように）グラフを追加するとか特徴があるようで。

http://d.hatena.ne.jp/foo22222/20110710/1309868771 も，ちゃんと書くとすれば以下のようにでも。

# ベータ関数の確率密度分布を描く.
# shape1=7, shape=15 の場合と shape1=1, shape2=1 の場合.
N <- 20
k <- 6
x <- seq(0, 1, length=10000) # ＠＠＠ これを追加
# "shape1=7, shape2=15"のベータ関数の確率密度分布.
plot(x, dbeta(x, shape1=k+1, shape2=N-k+1), ylim=c(0, 4), type="l") # ＠＠＠ 最初のグラフは plot で
#par(new=T) # ＠＠＠ こんなの使わない（きもちわるい）
# "shape1=1, shape2=1"のベータ関数の確率密度分布.
lines(x, dbeta(x, shape1=1, shape2=1), col="red") # ＠＠＠追加のグラフは lines や points で
# カンマの後には1個の空白を！！！

ダメ出し：12.3 いろいろな信頼区間

http://d.hatena.ne.jp/foo22222/20110710/1309865390

トリッキーなプログラミングで短くしてみた

x <- 1:10
s <- x*3
n <- x*10
ans <- array(unlist(mapply(binom.confint, s, n, prior.shape1=1, prior.shape2=1)), dim=c(11, 6, 10))[c(2, 3, 5, 9), 4:6,]
matplot(t(apply(ans, 3, "+")), type="l", col=c("black","red","blue","green"), lty=1)

ダメ出し：順列の逆関数

http://d.hatena.ne.jp/ryamada/20110622

これはちょっとすごい（ひどい）ですね。
定義された InvSeq は order 関数の結果と同じになります。

s が長いと，悲惨な結果になります。下の例では 1 万倍以上の速度比。

> s <- sample(10000) # sample(1:10000, 10000)  などとする必要はない

> system.time(invs <- InvSeq(s))
   user  system elapsed
 22.910   3.197  25.971

> system.time(invs2 <- order(s)) # 一瞬で終わる
   user  system elapsed
  0.002   0.000   0.002

> all.equal(invs, invs2)
[1] TRUE

> 22910/2
[1] 11455

ダメ出し：ぐるぐる回るために

http://d.hatena.ne.jp/ryamada/20111002　について

Npt <- 10000
n <- 3
X <- matrix(rnorm(Npt*n), Npt, n)
X <- X/sqrt(rowSums(X^2))
Y <- rowSums(X[,1:(n-1)]*X[,2:n])

は

Npt<-10000
n<-3
X<-matrix(rnorm(Npt*n),Npt,n)

X<-X/sqrt(apply(X^2,1,sum))
#apply(X^2,1,sum)

Y<-rep(0,Npt)

for(i in 1:(n-1)){
    Y<-Y+X[,i]*X[,i+1]
}

より，20 倍速い

ダメ出し：与えられた条件から、推定する～信頼区間？～

相変わらず for がお好き

http://d.hatena.ne.jp/ryamada/20110928

p<- seq(0,1,0.00001)
pp<- rep(0,length(p))
for(i in 1:length(p)){
    pp[i]<- dbinom(4,15,p[i])
}

より，35 倍速い。dbinom もベクトル化されている。使わないという手はない。

p<- seq(0,1,0.00001)
pp <- dbinom(4, 15, p)

ダメ出し：最小全域木とエッジの長さの和

http://d.hatena.ne.jp/ryamada/20110825

どうも，作者は for がお好きなようだが。

Npt<-100
k<-4
X<-matrix(0,Npt,k)
for(i in 2:Npt){
    X[i,]<-X[i-1,]+rnorm(k)
}

は，以下と同じことをしているはずだ。乱数発生もまとめて行う方が効率的。Npt=10000 にして計測すると，この部分だけだと 10 倍速いけど，それは全体の 5% 程でしかないので，トータルでは違いはないということになる。

Npt <- 100
k <- 4
x <- matrix(rnorm(Npt*k), Npt)
x <- matrix(apply(x, 2, cumsum), Npt)

ダメ出し：パーミュテーションテストのシャッフルの準備

http://d.hatena.ne.jp/ryamada/20110824

説明がよくわからないけど，プログラムを読む限りでは，行列の列を指定して，列ごとの要素をシャッフルするようなプログラムと言うこと。

sample 関数は，"与えられた要素をサンプリング（シャッフル）する" のだ。与える要素は 1:n でもよいし，x[1], ..., x[n] でもよい。しかし，1:n をシャッフルして，それに基づいて，x をシャッフルする必要などない。つまり，x <- c(4, 6, 2, 8) をシャッフルするなら，sample(x) だけでよい。よって，プログラムは簡単至極なものになる。

shuffle <- function(M, s=NULL)
{
    if (is.null(s)) {
        s < 2:ncol(M)
    }
    return(cbind(M[, -s], sapply(s, function(i) sample(M[, i]))))
}

実行例

> set.seed(888)
> n1 <- 3
> n2 <- 4
> (M <- matrix(sample(n1*n2), n1, n2))
     [,1] [,2] [,3] [,4]
[1,]    1    7    3    6
[2,]    4    9    2   11
[3,]   12   10    8    5
> shuffle(M, s=2:4)
     [,1] [,2] [,3] [,4]
[1,]    1    7    8    5
[2,]    4   10    2   11
[3,]   12    9    3    6

ダメ出し：カーネル密度分布推定の中身を覗く

http://d.hatena.ne.jp/ryamada/20110820 は，ベクトルプログラミングで 90 倍速くなる

300 ポイント作成では system.time が測定しきれないので，3000 ポイントに増やして測定する。

Npt <- 30000
x <- y <- rep(0, Npt)
r <- 10
t <- runif(Npt)*2*pi/1
r2 <- r*(sin(t*200)+10)+rnorm(Npt)*0.1
x <- r2*cos(t)
y <- r2*sin(t)
x <- cbind(x, y)

ダメ出し：正八面体タイル

正八角形タイルか?

http://d.hatena.ne.jp/ryamada/20110717 も以下のようになる

速度はほぼ2,3倍。なぜにして，for (i ... が必要なのか線を濃くするだけ?無駄も多いし。for は大敵。

func <- function(ox, oy, or) # 円を描く関数
{
    lines(ox+xc*or, oy+yc*or)
}
depth <- 3
k <- 8
a <- 1
x0<-0
y0<-0
as <- seq(0, 2*pi, length=100)
xc <- cos(as)
yc <- sin(as)
ts <- seq(0, 2*pi, length=8)*(k-1)/k
plot(0, 0, type="n", xlim=c(-1.5, 1.5), ylim=c(-1.5, 1.5))
b <- a/2*tan(pi/k)
c <- sqrt((a/2)^2+b^2)
r <- c*tan(2*pi/k)
R <- sqrt(r^2+c^2)+c
func(x0, y0, a)
sapply(ts, function(theta) func(R*cos(theta), R*sin(theta), r))

ダメ出し：双曲線幾何の直線その２

http://d.hatena.ne.jp/ryamada/20110713 のプログラムも読みにくく，効率が悪い。しかも途中で負数の平方根を求めるということになるのが不都合。

下のように書き直せば，100倍ほど速くなる。

func <- function(ox, oy, or) # 円を描く関数
{
     lines(ox+xc*or, oy+yc*or)
}
r1 <- 1
x2 <- 1.5
y2 <- 0
r2 <- 2
t <- seq(0, 2*pi, length=200)
xc <- cos(t)
yc <- sin(t)
plot(0, 0, type="n", xlim=c(-1.98, 3.5), ylim=c(-1.98, 3.5), asp=1)
func(x1, y1, r1)
func(x2, y2, r2)
k <- (x2^2+r1^2-r2^2)/(2*x2*xc)
r3 <- k^2-r1^2
mapply(function(x, y, z) if(z > 0) func(x*cos(y), x*sin(y), sqrt(z)), k, t, r3)

ダメ出し：双曲線幾何の直線

http://d.hatena.ne.jp/ryamada/20110709 のプログラム

もとのプログラムは，読みにくいこともさることながら，無駄なことをやっている。以下のように簡単になる。

## hyperbolic disk

t <- seq(0, 2*pi, length=100)
xc <- cos(t)
yc <- sin(t)
r <- tan(seq(from=0, to=pi/2, length=10))[-1]
func <- function(ox, oy, or) # 円を描く関数
{
    lines(ox+xc*or, oy+yc*or)
}
plot(0, 0, type="n", xlim=c(-1.5, 3), ylim=c(-1.5, 3), asp=1)
func(0, 0, 1)
sapply(c(r, -r), function(or) func(1, or, or))
abline(h=0)

短いプログラムはたいてい速い。案の定，元のプログラムの8倍ほど速くなっている。

Sweave-11

\include に対応するものはないけど \input に対しては，\SweaveInput がある。

\includeonly ができないのだけど，やむを得ない。

たくさん参加では，やはり統一性に欠ける？

> (mat <- matrix(1:24, 6))
     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    2    8   14   20
[3,]    3    9   15   21
[4,]    4   10   16   22
[5,]    5   11   17   23
[6,]    6   12   18   24

> subset(mat, mat[,1] %% 2)
 以下にエラー subset.matrix(mat, mat[, 1]%%2) :
   'subset' は論理値でなければなりません

> subset(mat, mat[,1] %% 2 == TRUE)
     [,1] [,2] [,3] [,4]
[1,]    1    7   13   19
[2,]    3    9   15   21
[3,]    5   11   17   23

mat[,1] %% 2　は 1 0 1 0 1 0 なので，他の場面では TRUE FALSE  TRUE FALSE  TRUE FALSE と解釈されることが多いのになあ？

いや，別に，ちゃんと記述すれば良いだけのことだから，文句は言わない。

 

Sweave-10

xtable と print.xtable の組み合わせで表を書くというのも，長ったらしいオプションであるとか，標準と違うデフォルトとかいうことで，これもまたラッパー関数を書くしかないかなと言うことで以下のようなものを作る。

ytable <- function(x, caption, label, align=NULL, digits=NULL, display=NULL,
                         table.placement="htbp", caption.placement="top",
                         math.style.negative=TRUE )
{
    library(xtable)
    xtable:::print.xtable(xtable(x, caption=caption, label=label, align=align, digits=digits, display=display),
                         table.placement=table.placement, caption.placement=caption.placement,
                         math.style.negative=math.style.negative)
}

.Rnw で

ytable(iris[1:10,], "iris data", "sss")

とか

x <- 1:10
y <- rnorm(10)
ytable(lm(y~x), caption="aaaaa", label="qqq")

caption と label は必須。もし，print.xtable の引数でデフォルトにしたいあるいはデフォルト値を変えたいものがあれば，この関数で定義しておく。

Sweave-9

図をフロートにする方法はSweave-3 に書いたが，\begin{figure} から始まる一連の指示を毎回書くのは面倒くさい。

表は xtable でそのままフロートになるのだから xtable のような関数を書いてやればよい。

xfigure <- function(fn, caption, label, location="htbp")
{
    cat(sprintf("\\begin{figure}[%s]\n", location))
    cat("\t\\centering\n")
    cat(sprintf("\t\\includegraphics{%s}\n", sub(".(pdf|PDF)", "", fn)))
    cat(sprintf("\t\\caption{%s}\n", caption))
    cat(sprintf("\t\\label{%s}\n", label))
    cat("\\end{figure}\n")
}

これを .Rnw で以下のように使う。

<<results=tex, echo=false>>=
fn <- "test.pdf"
pdf(fn, width=3, height=2, onefile=FALSE)
hist(rnorm(10000))
invisible(dev.off()) # invisible にするか，結果を何かに代入するかしないと余計な文字列が出力される
xfigure(fn, "ヒストグラム", "lab1")
@

毎回 invisible(dev.off()) と書くのが面倒なら，dev.off2 <- function invisible(dev.off()) とでもしておく。