裏 RjpWiki

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

プログラムを書く「お作法」

2013年01月28日 | ブログラミング

昔々,『プログラム書法』と『ソフトウェア作法』という本があったそうな。

異国の著者 Brian W.Kernighan と P.J.Plauger の著作を 木村 泉 氏が翻訳したものであったそうな。

「作法」は「作り方の法則」が本来の意図であろうが,木村氏はそれを「お作法(さほう)」つまり,「従わなければならない規則,あるいは人間の社会生活にかかわる多くの慣習やしきたりのうち,ふつうとくに起居動作,言語,身なりなどに関する正しいとされる方式を作法と呼ぶ」の観点でとらえたのがユニークといえたのだろうともいわれる。当時,外国語のコンピュータ用語を「日本語(漢字)」で表そう(アルゴリズムを算法,プログラムを算譜とか)というムーブメントのようなものもあったので,なんとなく心ひかれた人もいたそうな。

このブログで取り上げたたくさんの事例も,その観点に合致するものも多いかもしれない。言語環境,コンピュータ環境,コンピュータアーキテクチャから見れば,今となっては時代遅れ,見当違いの事例も多いかもしれないが,考慮すべき事例も数多く残っているかもしれない。それらを現代的な観点から振り返ることも意味のあることかもしれないなあ。なんちゃって。

コメント

そういうなら,そうしなさい

2013年01月28日 | ブログラミング

人柄ではなく『コード柄』あなたは大丈夫ですか?」に書いてあることだけど,

> この人は美しいコードを書くとか,行き届いたコードを書くとか。逆に,乱雑なコードを書く人もいれば,センスのないコードを書く人もいる…と。

> 「コード柄」は見る人が見れば一目瞭然。その人が「使える/使えない」を判断するうえで,職務経歴書を提出してもらうよりもずっと話が早いですよ,というエンジニアもいたのです

上から目線のお言葉ですが,そのサイトに掲載されているコードは,その基準を満たしていますか?

 

うわ。もっと,上から目線だわ。ごめんね。

コメント

アドバイス:問題文もレベルに応じて

2013年01月24日 | ブログラミング

http://gihyo.jp/dev/serial/01/codeiq/0007
第7回 t検定による問題解決,Rで実践できますか?~データサイエンティストの統計学─倉橋一成からの問題」に書かれていることなんだけど,

>    A群とB群⇒ 平均値が170cmの集団
>    C群⇒ 平均値が175cmの集団
>
> # A群のデータ生成
> set.seed(1)
> heightA <- 170 + 10*rnorm(100000)
> # B群のデータ生成
> set.seed(2)
> heightB <- 170 + 10*rnorm(100000)
> # C群のデータ生成
> set.seed(3)
> heightC <- 175 + 10*rnorm(100000)

まあ,それでもいいのだけど,rnorm には,mean と sd という引数があるんだから,

heightA <- rnorm(100000, 170, 10)

と書く方がスマートでしょう。

なお,アテーメー(当然)のことであるが,このようにして作成したデータセットの平均値と標準偏差は正確に 170 と 10 ではない。ちょうどにしたいなら,

heightA <- scale(rnorm(100000))*10+170

のようにすればよい。

なお,この準備に引き続いて,それぞれから 100 人をサンプリングせよという出題になっているけど,無限母集団から 100 人をサンプリングするのと,100000 人から 100 人をサンプリングすることの違いは意味がないと思うが??

問2のプログラムで for を使うテンプレートが書かれているのだけど,適切とは思われない。

「スキルを試すコードパズル」ということなので,出題様式が初心者レベルというのでは,なめられる恐れがある。

コメント

shiny でウェッブアプリ

2013年01月21日 | ブログラミング

二群の平均値の差の検定」に関するシミュレーションを行う shiny アプリ

http://glimmer.rstudio.com/rder/difference-of-two-means/ のように指定する(指定されたディレクトリに server.R, ui.R を置く)

曲線は母集団の密度関数

信頼区間は,母平均の 95% 信頼区間

母集団の密度関数はデータの分布を表している。

データの分布の重なりが相当大きくても,信頼区間が重なっていても,平均値の差についての検定結果は有意(二群の母平均には差がある)である。

「信頼区間が重ならないときに有意だ」というのはウソである。

そのほか,サンプルサイズ,平均値,標準偏差を様々に変化させて,P 値がどうなるかを観察し,考察せよ。

コメント

アドバイス:オンラインヘルプはよく読もう

2013年01月16日 | ブログラミング

Rでlmコマンドの結果からStd. Errorを取得する」に書いてあることについて

回帰分析の係数の区間推定を求めたいときがあり,そのためには,回帰係数の標準誤差を求める必要があるということですが。

その後,自分で計算式を書いて求めているのだけど,

回帰係数の信頼区間を求めるための confint という関数があるよということがオンラインヘルプに書いてある(See Also の項)。

> r <- lm(Sepal.Length~., iris[-5]) # expressionは任意の推定式です
> coef(summary(r))[, 2]
 (Intercept)  Sepal.Width Petal.Length  Petal.Width
  0.25077711   0.06664739   0.05671929   0.12754795
> # 例えば3番目の係数の95%信頼区間を求めたいときは、以下のようにしてください。
> a <- 0.05/2 # 両側検定
> df <- summary(r)$df[2] # 自由度
> b <- coef(r)[3] # 3番目の係数
> se <- coef(summary(r))[, 2][3] # 3番目の標準誤差
> sprintf("%.3f(95%%信頼区間%.3f~%.3f)", b, b-se*qt(a, df), b+se*qt(a, df))

[1] "0.709(95%信頼区間0.821~0.597)"

信頼区間を求めるには,

> confint(r)
                  2.5 %     97.5 %
(Intercept)   1.3603752  2.3516197
Sepal.Width   0.5191189  0.7825554
Petal.Length  0.5970350  0.8212289
Petal.Width  -0.8085615 -0.3044038

だけでよい。summary する必要も,その結果から SE を抽出する必要も,結果を計算する必要もない。返されるのは matrix

コメント

ダメ出し:何がやりたいのかを簡明直裁に書こう

2013年01月08日 | ブログラミング

MeanDecreaseGiniが大きい変数の箱ひげ図を描く」の記事についてだけど...

ベクトル要素の大きい順に 1, 2, ... と順序づけを行いたかったときに,「order 関数を二重に使えばよい」と twitter で教えてもらったということのようだ。

> x <- c(5,2,1,6,4,7,9,3,8)
> order(order(x, decreasing = T))
[1] 5 8 9 4 6 3 1 7 2

この問題をすなおに解こうとすると,「順位付け」というキーワードが浮かぶだろう。「ランク付け」ですね。そして,文字通り rank という関数があるんです。

? order をやると,See Also の所に, sort, rank, xtfm の 3 つが挙げてある。rank 関数は,同順位をどのように処理するかも指定できる。

で,やってみると(やってみなくてもわかるが)

> rank(x)
[1] 5 2 1 6 4 7 9 3 8

おや,結果が違う。それは,小さい順に順位付けされるから。では,大きい順に順位付けしたいためには?

ベクトル要素の符号が混じっていないなら(正の要素ばかりとか負の要素ばかりとか),マイナスをつける。

> rank(-x)
[1] 5 8 9 4 6 3 1 7 2

ほらできた。

しかしですね,ここで喜んではいけない。

筆者は,rk を求めて,最終的に n6 を得るためになにか不思議なことをやっている。

要素の値が大きい順にその要素の番号を取り出そうというのだから,そんな大変なわけがない。というか,order 関数の戻り値そのままじゃぁないの?

ということで,1 行で書けて,結果は同じだった。ということは,筆者は他人にわかりにくい,無駄なことをやったということか。

> rk <- order(order(rf.model$importance, decreasing = TRUE))
> n <- seq(1, nc2, 1)
> n2 <- cbind(rk, n)
> n3 <- data.frame(rk = n2[, 1], n = n2[, 2])
> n4 <- n3[order(n3$rk), ]
> n5 <- n4$n
> n6 <- n5[1:x]
> n6
[1] 52 53  7 16 55 56 21 25 57

> n6 <- order(rf.model$importance, decreasing = TRUE)[1:x]
> n6
[1] 52 53  7 16 55 56 21 25 57   同じだ

コメント

いろんなことをしてくれる関数が山ほどあるが...

2013年01月07日 | ブログラミング

reshape2 でデータ整形」に書いてあることについてなのだが

library(reshape2)
aq.m <- melt(airquality, id=c("Month", "Day"), na.rm=TRUE)
dcast(aq.m, Month ~ variable, mean)

によって

  Month    Ozone  Solar.R      Wind     Temp
1     5 23.61538 181.2963 11.622581 65.54839
2     6 29.44444 190.1667 10.266667 79.10000
3     7 59.11538 216.4839  8.941935 83.90323
4     8 59.96154 171.8571  8.793548 83.96774
5     9 31.44828 167.4333 10.180000 76.90000

更に,続いて,

library(plyr)
dcast(aq.m, Month~variable, mean, subset=.(variable=="Ozone"))

により

  Month    Ozone
1     5 23.61538
2     6 29.44444
3     7 59.11538
4     8 59.96154
5     9 31.44828

が得られると。

しかし考えても見よう,これくらいのことをするのに何故新たなパッケージを使わなければならないのか

aggregate 1 つで十分だ。長いって?複雑だって?そんなこたぁねぇなぁ~~~。

> aggregate(airquality[,1:4], by=airquality["Month"], mean, na.rm=TRUE)
  Month    Ozone  Solar.R      Wind     Temp
1     5 23.61538 181.2963 11.622581 65.54839
2     6 29.44444 190.1667 10.266667 79.10000
3     7 59.11538 216.4839  8.941935 83.90323
4     8 59.96154 171.8571  8.793548 83.96774
5     9 31.44828 167.4333 10.180000 76.90000
> aggregate(airquality["Ozone"], by=airquality["Month"], mean, na.rm=TRUE)
  Month    Ozone
1     5 23.61538

2     6 29.44444
3     7 59.11538
4     8 59.96154
5     9 31.44828

コメント (1)

ダメ出し:道具を作る場合には

2013年01月07日 | ブログラミング

全ての変数の分布を箱ひげ図にする」に書かれていることについて

提示されているプログラムは,作者本人も認めていること以外にもいろいろ不都合がある

道具は,とにかく「汎用的でなくてはならない」

そのためには,特定の状況を仮定したものではダメ。例えば,分析に使うデータフレーム名や,列数,列の構成を一般的なものではだめだというのは「いろはのい」。iris[,i] なんてのがでてきてはだめ。その対処は簡単で,関数にすればよいのだ。関数にするだけで,関数にするにはどうしなければならないかの対処をするだけで(それが難しいかも知れないが),自ずと汎用的な道具ができる。

そのほかにもいろいろあるが,あれこれ書くより,プログラム例を挙げる方が速い。どこがどう違うかは,みればわかる(ただし,真意が伝わらない部分もあるだろうが)。

boxplot2 <- function(d, g, nr=2, nc=2, width=500, height=400) {
    m <- ncol(d)
    png("fig%03i.png", width=width, height=height)
    layout(matrix(1:(nr*nc), nr, nc, byrow=TRUE))
    for (i in seq_len(m)) {
        boxplot(d[, i] ~ g, main = names(d)[i])
    }
    layout(1)
    dev.off()
}

boxplot2(iris[, 1:4], iris[, 5], nc=1)

boxplot2 関数の第 1 引数は,箱ひげ図を描く変数だけを含むデータフレーム。飛び飛びの列もちゃんと指定できる。

第 2 引数は,グループを表すベクトル。通常は,対象とするデータフレームの一つの列を指定する。

第 3,4 引数は,一つのファイルに何行何列の図を描くかを指定する。nr*nc が第 1 引数の列数の倍数である必要はない(最後のファイルには nr*nc 個未満の図が描かれるだけ)。

第 5,6 引数は,一つのファイルの画像サイズ(個々の箱ひげ図のサイズではない)。例では png を使ったが,pdf でもよい。それぞれ特徴があり,サイズの指定法も異なる場合がある。個人的には,pdf で onefile=TRUE がお勧めだ。

 

コメント