裏 RjpWiki

文字通り,RjpWiki の裏を行きます
R プログラム コンピュータ・サイエンス 統計学

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

2011年10月31日 | ブログラミング

どうも,この一派(?)は,プログラムはスペースを使わずギッチリ詰めて書くとか,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 いろいろな信頼区間

2011年10月31日 | ブログラミング

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)

コメント

ダメ出し:順列の逆関数

2011年10月24日 | ブログラミング

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

コメント

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

2011年10月24日 | ブログラミング

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 倍速い

コメント

ダメ出し:与えられた条件から、推定する~信頼区間?~

2011年10月24日 | ブログラミング

相変わらず 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)

コメント (1)

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

2011年10月24日 | ブログラミング

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)

コメント

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

2011年10月24日 | ブログラミング

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

コメント

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

2011年10月24日 | ブログラミング

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)

コメント

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

2011年10月24日 | ブログラミング

正八角形タイルか?

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))

コメント

ダメ出し:双曲線幾何の直線その2

2011年10月24日 | ブログラミング

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)

コメント

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

2011年10月24日 | ブログラミング

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

2011年10月17日 | RとLaTeX

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

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

コメント

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

2011年10月16日 | ブログラミング

> (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

2011年10月14日 | RとLaTeX

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

2011年10月14日 | RとLaTeX

図をフロートにする方法は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()) とでもしておく。

コメント