bdiag

無いものは作るという習慣がないのかな
あるかないかわからないものを探そう，教えてというのがなんだかなあ
以下のようなのでどうだ
わざわざわかりにくくしておいた(^_^;)
bdiag2 <- function(L)
{
a <- sapply(L, function(x) {x <- as.matrix(x); return(c(nrow(x), ncol(x)))})
l <- apply(a, 1, cumsum)
n <- nrow(l)
m <- array(0, dim=rowSums(a))
mapply(function(x, i, j) m[i, j] <<- x, L, mapply(seq, c(0, l[-n,1])+1, l[,1]), mapply(seq, c(0, l[-n,2])+1, l[,2]))
return(m)
}
L <- list(matrix(seq(4),2,2),matrix(seq(9),3,3),matrix(seq(25),5,5))
bdiag(L)
bdiag2(L)
mlist <- list(1, 2:3, diag(x=5:3), 27, cbind(1,3:6), 100:101)
bdiag(lapply(mlist, as.matrix))
bdiag2(lapply(mlist, as.matrix))

閉曲線

久保さん曰く
> 下の図のように不定形のセン (赤色) のことを敵国語で何と呼べばよいか?


久保さんが自答

> センだから line だろうと思うのはまちがいで …… line は「直線」「線分」のたぐいをあらわす (線分は正確には (line) segment だけど)． 曲線一般は curve なんだけど， 上のようなぐるっと一周まわって閉じてるモノに適用するのもなぁ …… まちがいではないけれど， いまいちなかんぢがする．

> ということで， 数学用語から借用すると 閉包 (closure) ということになる． しかし一般用語としての closure は閉鎖とかそういう意味だしなぁ …… と試行錯誤しつつ画像検索してみると， closed contour あたりが無難そうだ ( closed region てのもあるけど)．

======

私は今まで，closed curve だと思っていた
google でも用例は多い

Rcmdrでの複数変数の選択

> Rcmdrで重回帰分析やクラスター分析をする場合、変数を複数選択するところが出てきて、マウスでクリックすることで複数選択しますが、１つ置きに選択することができません。（連続したものしか選択できない）。一つ置きに選択する方法があるのでしょうか？

今更何を聞いているのか。他のソフト（Excel や Word!!）の場合だって，複数選択は ctrl キーを押しながらマウスクリックではないのか？ shift キーを押しながらマウスクリックというのは知っているようなんだが，知識が浅薄。

sd(list())は？

x がリストのとき，sd(x) はどうなるかということ

> sd
function (x, na.rm = FALSE)
{
if (is.matrix(x))
apply(x, 2, sd, na.rm = na.rm)
else if (is.vector(x))
sqrt(var(x, na.rm = na.rm))
else if (is.data.frame(x))
sapply(x, sd, na.rm = na.rm)
else sqrt(var(as.vector(x), na.rm = na.rm))
}


なので，sd(list()) はどうなるか。
事前にクラスを確かめておこう
> x <- list(d=1:10)
> class(x)
[1] "list"
matrix でも vector でも data.frame でもないので最後の選択肢
else sqrt(var(as.vector(x), na.rm = na.rm))
が実行されるが，これはエラーになる
> sd(x)
エラー： is.atomic(x) is not TRUE

> sqrt(var(as.vector(x), na.rm=TRUE))
エラー： is.atomic(x) is not TRUE

> var(as.vector(x), na.rm=TRUE)
エラー： is.atomic(x) is not TRUE

is.atomic がなんで TRUE じゃないのかと見てみると

> as.vector(x)
$d
[1] 1 2 3 4 5 6 7 8 9 10

> class(as.vector(x))
[1] "list"

そりゃ，list ですよね
ということで，
> sd(as.vector(x)$d, na.rm=TRUE)
[1] 3.027650
> var(as.vector(x)$d, na.rm=TRUE)
[1] 9.166667
>

sd の実装が十分なのかどうかはコメントしづらい

二項分布の上側確率

母比率 0.5，試行回数 10 の二項分布で，x が 9 以上になる確率を求めるのに，

> pbinom(9, 10, 0.5, lower.tail=FALSE)
[1] 0.0009765625

としてはだめだ。

> pbinom(8, 10, 0.5, lower.tail=FALSE)
[1] 0.01074219

でなくてはならない。lower.tail=TRUE と，以下を考え合わせれば分かるが，うっかりすると最初に挙げたようにしてしまいそう。

> cbind(0:10, x, cumsum(x))
x
[1,] 0 0.0009765625 0.0009765625
[2,] 1 0.0097656250 0.0107421875
[3,] 2 0.0439453125 0.0546875000
[4,] 3 0.1171875000 0.1718750000
[5,] 4 0.2050781250 0.3769531250
[6,] 5 0.2460937500 0.6230468750
[7,] 6 0.2050781250 0.8281250000
[8,] 7 0.1171875000 0.9453125000
[9,] 8 0.0439453125 0.9892578125
[10,] 9 0.0097656250 0.9990234375
[11,] 10 0.0009765625 1.0000000000

R 2,7 の新機能など

いろいろ変更点が多い。
choose の定義なども，びっくりする。help を見るべし。
> choose(-11,4)
[1] 1001
> choose(-11,3)
[1] -286
> ? choose
> choose(12.3, 5) # この例なんか，びっくり^2 だな
[1] 920.8748

グラフィックファイルのデフォルト

2.7 から，ps ではなく pdf が標準になった。pdf 関数の file 以外の引数は pdf.options で設定できる。dev.copy2pdf とかいろいろ ps 関連の関数に対して pdf 関連の関数が用意された。

紅一点

ほとんど数値であるデータファイルからデータを読んでいて，数値じゃないものが出てきたら，エラーにしなくても良いけど，警告くらい出したらどうだろうかね。
----- file start
a, b
8, l
9, 2
1o, 3
11, 4
----- file end
というデータファイルを，データフレームに読み込んで，a, b の要素を調べてご覧あれ。

二次データのt検定

R の t.test は，元データがないと使えないと思っている人がいる。でも，裏技（？）がある。
ヒントは，指定されたサンプルサイズで，指定された平均値と標準偏差を持つデータが作れるということ。t検定は本質的にはサンプルサイズと平均値と標準偏差しか使わない。個々のデータがどういう値を取るかとか分布は正規分布かなどということを問わないのだ。
ということで，サンプルサイズが n1, n2，平均値がm1, m2，標準偏差がs1, s2 ということしか分かっていないときに t.test を使うには，，，
t.test(scale(rnorm(n1))*s1+m1, scale(rnorm(n2))*s2+m2)　これだけでよいのだ～～っ。

> t.test(scale(rnorm(31))*3.212+154.546, scale(rnorm(16))*4.387+153.698)

Welch Two Sample t-test

data: scale(rnorm(31)) * 3.212 + 154.546 and scale(rnorm(16)) * 4.387 + 153.698
t = 0.6843, df = 23.547, p-value = 0.5005
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-1.712228 3.408228
sample estimates:
mean of x mean of y
154.546 153.698

念のために書いておくと，何もわざわざ架空のデータベクトルなど作らなくても，各群のサンプルサイズ，平均値，標準偏差の，合計6つの引数からt検定を行う関数を作ればよいのだよ。

*apply関数一族

tapplyを使ったグラフ作成 (2008-01-11 (金) 11:47:06) の結論は

AGE<-c(30, 80, 90, 77, 65, 33)
TYPE<-factor(rep(c("male", "female"), each=3), levels=c("male", "female"))

myhist <- function(x) {
hist(x, main=paste("Histogram of TYPE=", levels(TYPE)[substitute(x)[[3]]]))
}

tapply(AGE, TYPE, myhist)

ということですね。

関数の中で大域変数を参照するのは避けた方がよいと思われ，以下のように書いてみる。。。

myhist <- function(x, g) {
hist(x, main=paste("Histogram of TYPE=", g))
}

a <- mapply(myhist, split(AGE,TYPE), levels(TYPE))

他のデータを処理するときも，関数の中を弄らなくてもよいというメリットあり

b <- mapply(myhist, split(iris[,1], iris[,5]), levels(iris[,5]))

引き数名の省略

「引き数名は他の引数と識別できるかぎり省略可能」と説明されてきたが，例外もある。
「... の後に出てくる引数は完全に一致しなければならない」と。
x <- rnorm(10)
cut(x, break=c(-10, -1, 0, 1, 10))
は，シンタックスエラーになる。
breaks でないといけない。
あるいは，引き数名を省略する。
cut(x, c(-10, -1, 0, 1, 10))
===========================
> x <- rnorm(10)
> cut(x, break=c(-10, -1, 0, 1, 10))
エラー： syntax error
> cut(x, breaks=c(-10, -1, 0, 1, 10))
[1] (-1,0] (0,1] (0,1] (-1,0] (0,1] (0,1]
[7] (0,1] (0,1] (0,1] (-10,-1]
Levels: (-10,-1] (-1,0] (0,1] (1,10]
> cut(x, c(-10, -1, 0, 1, 10))
[1] (-1,0] (0,1] (0,1] (-1,0] (0,1] (0,1]
[7] (0,1] (0,1] (0,1] (-10,-1]
Levels: (-10,-1] (-1,0] (0,1] (1,10]

裏 RjpWiki

> WindowsXPでR2.6.1を使って一般化線形モデルによる解析を行っています。
> データを読み込んでさぁ解析といったところで、次のようなメッセージがでました。

> x=read.table("a.txt",header=T)

> result=glm(A~.,data=x)
> エラー： サイズ 294.8 Mb のベクトルを割り当てることができません

その後，以下もやったが

> memory.size(T)
> [1] 801.5625

なぜだか教えてというのだが。

memory.size のヘルプくらい読めばよいのに。
そうすれば，memory.size(FALSE) もやってみれば，現在どれくらいメモリを使っていて，あと294.8MB 使おうとしてだめだった理由もわかるだろうに（1GB積んでいるというのだから，引き算なり足し算なりすればわかるだろう）

−−−−−コメント

以下でわかるだろうと言うこと
> memory.size(TRUE)
[1] 13.1875
> memory.size(FALSE)
[1] 10.81405
> memory.limit()
[1] 246.4844
> gc()
used (Mb) gc trigger (Mb) max used (Mb) Ncells
138623 3.8 350000 9.4 350000 9.4 Vcells 81410 0.7 786432 6.0 472989 3.7

久保さんのデータ

http://hosho.ees.hokudai.ac.jp/~kubo/log/ の2008年1月6日のデータであるが。
元のデータに三次曲線を当てはめて，変化率の変化率が正から負へ向かうポイントを推定する。得られた三次曲線を二度微分してやって解を求めるとめでたく20 近辺の値になる。
というのは，久保さんが作ったデータがこの3字曲線に基づいているからそうなったんだろうなあと思う。
なお，図から座標を読み取るプログラムを用いた。この場合はデータが整数値なので，もとのデータは正確に再現できているだろう。
df <- structure(list(x = 1:50, y = c(165, 203, 200, 180, 153, 148, 181, 152, 173, 160, 154, 130, 171, 156, 162, 149, 163, 184, 151, 187, 187, 170, 165, 181, 193, 213, 207, 208, 204, 191, 217, 220, 202, 208, 182, 187, 176, 178, 145, 154, 136, 132, 115, 114, 98, 79, 78, 67, 60, 46)), .Names = c("x", "y"), class ="data.frame", row.names = c(NA, -50L))
ans <- lm(y ~ x+I(x^2)+I(x^3), df)
plot(df)
lines(df$x, ans$fitted.values, col="red")
summary(ans)
par <- (ans$coefficients*0:3)[-1]
par <- par*0:2
-par2[2]/par2[3]

encoding によらない入力

> データファイルのエンコーディング（UTF-8とSHIFT_JISなど）を区別せずに入力したい（若干修正）

そんなに色々なエンコーディングのファイルを入力する必要はないと思われるので，さして重要な要求ではないとは思うが，プログラム開発者の立場からは，そのような方法があればありがたいな。

−−−−−コメント

エディタなんかは，テキストのエンコーディングを判定してから開くのが当たり前なんだから，Rだって，そのような仕様にするのは簡単なことなはずなんだけ ど（私にはできないが）

na.print

print メソッドの引数 na.print の挙動がわからない。

data.frame の中の NA は na.print をどうしようとも NA としか表示されない。
vector や matrix の場合には期待される通りなんだが。ヘルプを見ると print.data.frame には na.print という引数はないようでもあるのだが。
> v <- c(1,4,NA,6,9)
> print(v, na.print="--")
[1] 1 4 -- 6 9
> m <- matrix(v, 5)
> print(m, na.print="--")
[,1]
[1,] 1
[2,] 4
[3,] --
[4,] 6
[5,] 9
> df <- as.data.frame(m)
> print(df, na.print="--")
V1
1 1
2 4
3 NA
4 6
5 9