裏 RjpWiki

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

応用力の養成が必要

2011年03月29日 | 統計学

「一を聞いて十を知る」という言葉がある。でも,それは簡単にはできないことだろう。一を聞いて,過去の知識・経験と照らし合わせて,さらにそこからスタートして色々調べて・試行錯誤して,「ああ!そうなのか!」,「おお!そうなるのか!」という経験をするのが,「学問」であり「修行」なのだろう。教えられたことだけを理解して(理解したつもりになって),教えられたことだけができるようになっても,それは「人間ではない,ロボットやコンピュータに過ぎない」。

ネットで色々質問する人がいて,それに対して答える人がいて,その関係を覧ているに付け,「この質問をした人は,一を聞いて十を知る人か否か」というのが,何となく分かる気がする。

コメント

「順列検定」って何?

2011年03月29日 | 統計学

あるブログで,近々刊行される「Rによる計算機統計学」の翻訳がでたらめだ(?)ということが縷々書かれていたが,その続報で,出版社のホームページを見てみると,目次の中に見慣れない用語があった。

「順列検定」というのが載っているのだけど,原語は permutation test。確かに,permutation は 順列だけど,permutation test は普通は「並べ替え検定」とされているのではないかな?

ググってみると,順列検定は 190 件,並べ替え検定は 2500 件。多ければ良いというわけではないが,大勢に従うのが吉かな。

そのほか,「確率変数を生み出す方法」という章があるが,これは,「確率変数を生成する方法」かな?

キャッチコピーに,「解析についてはRにより、統計学的コンピューティングの事例によるアプローチで計算統計の古典的な主要問題を解説した。」とあるが,「計算統計」はおかしいだろう。「計算機統計学」かあるいは「統計計算」ではないか?(たぶん前者だろうが)

大学院生かなにかを動員して,翻訳したのかなあと疑う。

コメント

ブートストラップとは

2011年03月04日 | ブログラミング

「ブーツストラップ」と誤解していたり,「ストラップ」に捕らわれて「編み上げ靴の紐」などと誤解している向きが多いが,本当は,靴(特に編み上げ靴?)のかかと部分にある「つまみ皮」のこと。

コンピュータが動くにはプログラムが必要だが,昔のコンピュータは電源を入れた直後にはメモリには何のプログラムも入っていない。そこで,スイッチの組合せなどで,「プログラムを読み込むためのプログラム」を設定する。当然短いものだから,それだけで十分な機能を持てるわけがない。そこで,少しましな機能を持つプログラムを読み込み,さらにそのプログラムはもっとましな機能を持つプログラムを読み込む。というようなことを何段階か重ねてすごいコンピュータプログラムをインストールする。この過程を「ブートストラップ」と呼ぶ。

なぜ,これをブートストラップと呼ぶかと言えば,「靴のブートストラップを思いっきり引っ張り上げれば,我が身は中に浮かぶ」という(元々はほら吹き男爵の話に出てくるというのだが,ここに書いてあるという証拠を出せる人がいないようなのだ)馬鹿馬鹿しいお話に端を発する。つまり,「コンピュータが自分を動かすためのプログラムを自分で生成する」みたいな点が似ていると言うこと。

コメント

R プログラムの抽出

2011年03月03日 | ブログラミング

ホームページなどで R プログラムが紹介されているとき,プログラムとその結果が混在している R コンソールをコピー&ペーストされたものが載っていることがある。

実際に試してみるとき,キー入力すれば済むことだけど面倒だし,プログラムの先頭の > や + を除いてコンソールにペーストするのもかなり面倒。

そんなときのために以下のような関数を作っておく。

使い方は,

(1) 関数を定義(source などで読み込む。環境を保存しておけば,毎回読み込む必要はない)

(2) Web ページなどの該当部分をコピー

(3) e() を実行

(4) R コンソールにペースト

注:ときどき,最終行がコピー&ペーストされないことがある

e <- function()
{
    con <- pipe("pbpaste", "r")
    txt <- readLines(con)
    close(con)
    con <- pipe("pbcopy", "w")
    for (str in txt) {
        if (grepl("^> |^\\+ ", str)) {
            writeLines(sub("^..", "", str), con)
        }
    }
    close(con)
}


コメント

ヒストグラムも汚い

2011年03月03日 | ブログラミング

http://scienceoss.com/2008/03/ に ggplot2 でヒストグラムを描く例がある。

library(ggplot2)
ggplot(iris)+aes(x=Sepal.Width)+geom_histogram()+facet_grid(Species~.)

というように,簡単な指定で描けるが,改めてみるまでもなく汚い。

灰色の地に灰色の棒はみじめ。プログラムを実際に動かしてみると,私の環境では,棒は真っ黒。
まあ,棒の色や地の色は変更することもできるのでどうということはないのだけど。

困るのは,ヒストグラムの階級が細かすぎるのと,きれいな数値でないこと。
まあ,自分で変えることはできるよというメッセージも出るけど,デフォルトがこれではちょっと困る。
このデータの最小値は2なのだけど,2のデータがどこに出るかと言えば「~2」のところにある。
欧米では日本と違って,「○○より大きく,□□以下」という階級定義なので,まあ,やむを得ない。
これも,階級の終わりがオープンかクローズかの指定もできるだろうけど,もはや,やる気が失せる。

さて,以下のようなプログラムで描いてみるとどうかな。

layout(matrix(1:3), 3)
old <- par(mar=c(3, 3, 0.1, 1), mgp=c(1.8, 0.6, 0))
sw <- split(iris$Sepal.Width, iris$Species)
xlim <- range(pretty(iris$Sepal.Width))
breaks <- seq(xlim[1], xlim[2], by=diff(xlim)/20)
ylim <- c(0, max(sapply(sw, function(x) table(cut(x, breaks, right=FALSE)))))
for (i in names(sw)) {
    a <- hist(sw[[i]], breaks=breaks, xlim=xlim, ylim=ylim, right=FALSE,
                col="aliceblue", main="", xlab=paste("Sepal.Width", i, sep=":"))

コメント

2 のべき乗の先頭桁

2011年03月02日 | ブログラミング

2のべき乗数,2^0, 2^1, 2^2, 2^3, ... の先頭桁を考える。

一見,1, 2, 4, 8, 1, 3, 6, 1, 2, 5, 1, 2, 4, ... と長さ10の周期で循環しそうに見えるがそれは違う。

先頭桁に,7,9 は出てこないように思うとそれは違う。

以下のプログラムで出力される行方向優先の行列が,先頭桁が何であるかを示している。

library(gmp)
func <- function(x)
{
    return(as.integer(unlist(strsplit(as.character(x), ""))[1]))
}
x <- as.bigz(1)
n <- 140
f <- integer(n)
for (i in 1:n) {
    x <- x*2
    f[i] <- func(x)
}
matrix(f, ncol=10, byrow=TRUE)

結果

      [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
 [1,]    2    4    8    1    3    6    1    2    5     1
 [2,]    2    4    8    1    3    6    1    2    5     1
 [3,]    2    4    8    1    3    6    1    2    5     1
 [4,]    2    4    8    1    3    6    1    2    5     1
 [5,]    2    4    8    1    3    7    1    2    5     1
 [6,]    2    4    9    1    3    7    1    2    5     1
 [7,]    2    4    9    1    3    7    1    2    5     1
 [8,]    2    4    9    1    3    7    1    3    6     1
 [9,]    2    4    9    1    3    7    1    3    6     1
[10,]    2    4    9    1    3    7    1    3    6     1

最初に先頭桁が 7 になるのは 2^46,9 になるのは 2^53 である

> as.bigz(2)^45
[1] "35184372088832"
> as.bigz(2)^46
[1] "70368744177664"

> as.bigz(2)^52
[1] "4503599627370496"
> as.bigz(2)^53
[1] "9007199254740992"

コメント

不適切なグラフ in ggplot2 その2

2011年03月02日 | ブログラミング

以下のようなグラフを描く関数を書いてみるが,仕様がダメダメ(ダメダメグラフの仕様がまともなわけがないので)なので,改造する気にもならない。

プログラムは以下のようなもの

hedo <- function(x, labels=NULL, col=NULL, acc=200)
{
    count <- table(x)
    pos <- count/max(count)*0.9+0.1
    n <- length(pos)
    theta <- seq(pi/2, -pi*3/2, length=acc)
    x <- cos(theta)
    y <- sin(theta)
    theta2 <- seq(pi/2, -pi*3/2, length=n+1)
    plot(x, y, type="n", col="gray80", asp=1, axes=FALSE, xlab="", ylab="")
    for (fraction in seq(0.1, 1.0, length=11)) {
        lines(x*fraction, y*fraction, col="gray80")
    }
    text(-strwidth("0")*0.7, seq(0.1, 1.0, length=11)+strheight("0")*0.3,
         sprintf("%3.1f", seq(0, 1, by=0.1 )), pos=4, cex=0.7)
    if (is.null(col)) col <- heat.colors(n, alpha=0.6)
    for (i in 1:n) {
        theta <- seq(theta2[i], theta2[i+1], length=acc)
        theta <- c(theta, rev(theta), theta2[i])
        fraction <- c(rep(pos[i], acc), rep(0.1, acc), pos[i])
        lines(c(0.1, 1.0)*cos(theta2[i]), c(0.1, 1.0)*sin(theta2[i]), col="gray80")
        polygon(cos(theta)*fraction, sin(theta)*fraction, col=col[i])
        mid <- (theta2[i]+theta2[i+1])/2
        text(cos(mid)*1.1, sin(mid)*1.1, paste(labels[i], count[i], sep=":"))
    }
}
#hedo(sample(10, 100, replace=TRUE))
hedo(x, labels=LETTERS[1:10])

 

コメント

不適切なグラフ in ggplot2

2011年03月02日 | ブログラミング

「以下に示すような図を『Excel で描くにはどうすればよいか』をビデオで説明しているサイト http://www.excelcharts.com/blog/the-consultants-chart-revisited/ があるが,大変手間がかかるようだけど,ggplot2 でやれば簡単だよ」というサイト http://learnr.wordpress.com/2010/08/16/consultants-chart-in-ggplot2/ があるけど,そもそもこんな図は不適切。ヨイコハマネシテハイケマセン。

これは ggplot2 で描いたものだというが,セクターの面積がデータの大きさを反映していない。

本当は,以下のような棒グラフを描くべきなのだ。

しかしまあ,なんですね~。なぜこうも,派手な色使いをするんだろうか。

色にも,膨張色と縮小色があるんだから,その影響もあるだろうし,何よりも,眼がチカチカする。

それにまた,右端の凡例なんて,意味ねー!

まあ,必要ないけど ggplot2 を使ったプログラム。

set.seed(9876)
DF <- data.frame(variable = 1:10, value = sample(10,
    replace = TRUE))
DF
library(ggplot2)
ggplot(DF, aes(factor(variable), value, fill = factor(variable))) + geom_bar(width = 1)
last_plot() + scale_y_continuous(breaks = 0:10) +
    coord_polar() + labs(x = "", y = "") + opts(legend.position = "none",
    axis.text.x = theme_blank(), axis.text.y = theme_blank(),
axis.ticks = theme_blank())

コメント