「楽しい 統計学セミナー - 証拠を科学する方法 -」の「第3講」 において
サイコロを6回振って,6の目が出る「確率分布図」を描いてみる.
:
X軸のラベルを確率変数Xと合わせることができない・・・ といっているが...
この項は,前の記事に既に解があるが...
x <- numeric(7)
for (i in 0:6) {
x[i +1 ] <- choose(6, i) * (1 / 6)^(i) * (5 / 6)^(6 - i)
}
plot(x, type="l", xaxt ="n", xlab="X", ylab="probability")
axis(side=1, at=0:6)
X 軸のラベルを確率変数 X と合わせるには,axis(side=1, at=1:7, labels=0:6) とすればよいだけ。
ここでも,実際の x の添え字は 0~6 だが,R では 1 ~ 7 になっているという食い違いを,ちゃんと吸収できていないことが原因(前の記事での指摘通り)。
度数分布を折れ線で描く(度数多角形という)ということもあるが,離散変数の分布は barplot 関数で描く方がよい。
R のbarplot は,横軸のラベルを描いてくれないので自分で描く。
確率分布関数として dbinom を使うということも,以前指摘したとおり。
barplot(dbinom(0:6, 6, 1/6), ylab="Probability", xlab="x")
axis(1, 0:6*1.2+0.6, 0:6)
「楽しい 統計学セミナー - 証拠を科学する方法 -」の「第2講」 において
試行回数 N=10~100 に対する二項分布の確率関数を
x <- numeric(91)
for (i in 10:100) {
x[i - 9] <- choose(i, 3) * (1/10)^(3) * (1- (1/10))^(i - 3)
}
plot(x)
abline(v = 20.5, col="blue")
のように,書いている。10 ~ 100 を添え字として 1 ~ 91 を使おうということで numeric(91) とか for (i in 10:100) とか x[i - 9] とか,様々な注意が必要。しかし,このようなことはちょっと油断するとバグの混入原因。縦軸(確率)を納める変数名が x なのも混乱を増幅するが。
実際,掲載されたプログラムが
for (i in 10:91) {
x[i - 9] <- > for (i in 10:100) {
となっていたり,確率が最大になる index の値を図に示すと,index=20, 21 のときに最大になるというように見えるが,求める数値は 29 と 30 なのだ,のように,ぼろが出る。
添え字は,0 やマイナスの値にならない限りそのまま使う方がよい。
y <- numeric(100)
n <- 10:100
for (i in n) {
y[i] <- choose(i, 3) * (1/10)^(3) * (1- (1/10))^(i - 3)
}
plot(n, y[n])
abline(v = 29.5, col="blue")
なお,二項分布の確率関数は dbinom で計算できる
> x <- 10:100
> y <- dbinom(3, x, 1/10)
> plot(x, y, type="l", ylim=c(0, 0.25))
> (index <- which.max(y))
[1] 20
> abline(v = x[index] + 0.5, col="blue")
> (max.value <- y[index])
[1] 0.2360879322323426
「楽しい 統計学セミナー - 証拠を科学する方法 -」の「第1講」 において,
当初目指していた「(問題)サイコロを6回投げたとき,6の目が1回だけ出る確率は? 」であるが,これは,生起確率が 1/6 である事象が 6 回の施行中に 1 回だけ起きるということで,それは二項分布だ。
Prob = dbinom(1, 6, 1/6) = 0.4018776 で,これを確かめるシミュレーションプログラムは,前のプログラムをちょっと変更すると得られる。
なお,水平線を引くのは abline(h = value, ...) である。
system.time({
set.seed(12345)
n <- 10000
x <- matrix(sample(0:1, 6*n, replace=TRUE, prob=c(5/6, 1/6)), 6)
y <- cumsum(colSums(x)==1)/1:n
plot(y, type="l")
abline(h = dbinom(1, 6, 1/6), col="red")
})
barplot(dbinom(0:6, 6, 1/6), ylab="Probability", xlab="x")
axis(1, 0:6*1.2+0.6, 0:6)
「楽しい 統計学セミナー - 証拠を科学する方法 -」の「第1講」 において,シミュレーションプログラムが掲示されているが
system.time({
set.seed(12345) # 乱数の初期値を設定しておこう
n <- 10000 # プログラムの中で何回も出てくる数値は変数に代入しておいて使おう(一箇所で変更できる)
# for (i in 1:n) { # 等差数列を作るというのだが...
# x[i] <- 6*i
# }
x <- 1:n*6 # ベクトル演算はこんなときにも使える (1:n)*6 と書くかは,好みの問題
y <- numeric(n) # 1~n回の試行を格納するための変数を確保.
for (i in 1:n) {
a <-sample(x=1:6, size=x[i], replace=TRUE, prob=c(rep(1/6, 6)))
# y[i] <- (sum(a[a == 1])) / x[i] # 足してデータの個数で割るというのは...
y[i] <- mean(a == 1) # 平均値を求めるということ
} # ちなみに,(a[a == 1]) は冗長
plot(y, type="l")
abline(a = 1 / 6, b = 0, col="red")
})
このプログラムでは,「6 回サイコロを振る」ということを,毎回 1~x[i] 繰り返す。
もう一つの考え方は,1 ~ n 回繰り返し,その途中経過を見るということでもよい。もっとも,こちらの考え方では k 回の時点での結果は k-1 回までの結果に従属であるということではあるが(前述の方法では,各回の結果は独立)。
もう一つのポイントは,1の目が出るということに注目するのであるから,2~6の目が出ることをシミュレーションする必要はないということ。つまり 1 とそれ以外の目が 1/6, 5/6 で出るというシミュレーションをすればよいということ。TRUE/FALSE や,成功回数を数えるときに sum が直接使えるということを考えれば,1 の目がでることを 1,それ以外の目が出ることを 0 で表すとよい。つまり n 回の試行結果は x <- sample(0:1, n, replace=TRUE, prob=c(5/6, 1/6)) で得られ,k 回目までの成功数は cumsum(x) で得られ,成功率は cumsum(x)/1:n で得られることになる。ということで,纏めると
system.time({
set.seed(12345)
n <- 10000
x <- matrix(sample(0:1, 6*n, replace=TRUE, prob=c(5/6, 1/6)), 6)
y <- cumsum(colSums(x))/6/1:n
plot(y, type="l")
abline(a = 1 / 6, b = 0, col="red")
})
前のプログラムの実行時間は
ユーザ システム 経過
15.418 2.327 17.636
後ろのプログラムの実行時間は
ユーザ システム 経過
0.136 0.004 0.155
"Performing multiple t-tests on different variables between the same two groups" で,ラッパーを紹介しているが
余り頻繁に使わないラッパーだと,ラッパーの使い方自体忘れたりして,そのたびに思い出すためにソースを見たり(メモを見たり)ということになる。
たとえば,t.test の場合なども,必要になったらすぐに簡単に書けるようにしておくと吉。
kidney[c("time", "age", "frail")] の部分は,分析対象とする変数だけを含むデータフレームを作っているだけなので,たくさんの変数を指定する場合などは変数名を書くだけで疲れるので, kidney[c(2, 4, 7)] のような指定法による方が便利。この程度のプログラムなら,ラッパーをわざわざ作っておく必要もないだろう。
関数部分が複雑(たとえば,これらの結果をまとめて一つの表にし,LaTeX のソースを作るなど)になると,新たに関数を作る意味もあるだろう。もっとも,そのような関数はもはやラッパーではないけど。
> library(survival)
> data(kidney)
> sapply(kidney[c("time", "age", "frail")], function(x) t.test(x~kidney$sex))
time age
statistic -1.733522 -0.06756775
parameter 34.53963 32.03715
p.value 0.09192286 0.9465497
conf.int Numeric,2 Numeric,2
estimate Numeric,2 Numeric,2
null.value 0 0
alternative "two.sided" "two.sided"
method "Welch Two Sample t-test" "Welch Two Sample t-test"
data.name "x by kidney$sex" "x by kidney$sex"
frail
statistic 0.8431417
parameter 27.62198
p.value 0.4063915
conf.int Numeric,2
estimate Numeric,2
null.value 0
alternative "two.sided"
method "Welch Two Sample t-test"
data.name "x by kidney$sex"
> lapply(kidney[c("time", "age", "frail")], function(x) t.test(x~kidney$sex))
$time
Welch Two Sample t-test
data: x by kidney$sex
t = -1.7335, df = 34.54, p-value = 0.09192
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-124.761114 9.861114
sample estimates:
mean in group 1 mean in group 2
59.30 116.75
$age
Welch Two Sample t-test
data: x by kidney$sex
t = -0.0676, df = 32.037, p-value = 0.9465
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-8.342454 7.806740
sample estimates:
mean in group 1 mean in group 2
43.50000 43.76786
$frail
Welch Two Sample t-test
data: x by kidney$sex
t = 0.8431, df = 27.622, p-value = 0.4064
alternative hypothesis: true difference in means is not equal to 0
95 percent confidence interval:
-0.2442908 0.5857194
sample estimates:
mean in group 1 mean in group 2
1.310000 1.139286
ggplot2 の geom_bar の作図例で "Stacked bar charts" の凡例とグラフの順序が逆。色使いも下品だが。
ggplot(diamonds, aes(clarity, fill=cut)) + geom_bar()
普通に barplot で描いてみよう。ちゃんと順序がそろう。
ans <- xtabs(~cut+clarity, diamonds)
col <- c("#00a0e960", "#e4007f60", "#00994460", "#f3980060", "#0068b760")
old <- par(las=1, mar=c(6, 6, 1, 1))
barplot(ans, legend=rownames(ans), col=col, xlab="clarity", ylab="count")
par(old)
# いびつなサイコロを1000回投げて,出る目の傾向を調べてみる.
x <- numeric(1000)
for (i in 1:1000){
x[i] <- sample(x=1:6, size=1, replace=FALSE, prob=c(1,2,3,3,2,1)) # 3と4の目が出やすいサイコロ.
}
これは,以下のように一行で書く(余分なことはしない。for など使うと遅くなる)
x <- sample(x=1:6, size=1000, replace=TRUE, prob=c(1,2,3,3,2,1))
> いびつなサイコロを1000回投げて,出る目の傾向を調べてみる
:
> あれれ・・・,出た目1と2のバーの間に隙間がない・・・?
その前に,1の棒は1の目盛りの右にあり,2の棒は2の目盛りの左にあるというおかしな状況はどう思います?なぜそのような図になるのかはオンラインヘルプに書いてある(かな?)。
right
logical; if TRUE, the histogram cells are right-closed (left open) intervals.
If right = TRUE (default), the histogram cells are intervals of the form (a, b], i.e., they include their right-hand endpoint, but not their left one, with the exception of the first cell when include.lowest is TRUE.
For right = FALSE, the intervals are of the form [a, b), and include.lowest means ‘include highest’
途中をすっ飛ばして結論を述べれば,hist は連続変数の分布図を描くためにあるのであって,離散変数は barplot で描くべきなのです,ということ。
barplot(table(x))
どうしても hist で描きたいというなら,
hist(x, breaks=1:7-0.5)
とすれば,見た目はもっともらしいものができるが...
「ggplot2 bar_plot with label」で,
「エクセル脳の皆さんに贈る」などといっているが,あなたこそ「ggplot脳」なのでは?
提示されたグラフは統計学的には不適切きわまりない。不適切なグラフに,どんなラベルを付けようが,不適切なグラフが適切なグラフになるわけはない。
適切なグラフは box-plot である。それが,「エクセル脳の皆さん」にも,「ggplot脳の皆さん」にも受け入れがたくても,仕方ない。
左のグラフを描くには,以下のようにする。
実に単純だ。なお,言うまでもないが,左の図の太い水平線は,平均値ではなく中央値だ。
boxplot(Sepal.Length~Species, data=iris)
median <- by(iris$Sepal.Length, iris$Species, median)
text(1:3+0.35, median, median, pos=4, xpd=TRUE)
「TokyoR25」にて,library(compiler) の威力を見せつけようと??しているようだが,compile したときとしないときの比較ではなくて,単に標準誤差はサンプルサイズの平方根に反比例するという統計学上の当たり前の話をなぞっているだけのような。
ここでも,もう何回となく取り上げてきて,うんざりということだが,プログラムを書く一人一人についてはうんざりでもなんでもなく,え~そんなことがあるの?ということなんだろうねえ。
元のプログラムは,シミュレーションごとに標本を生成している。それは,ダメだっっていってるでしょ。まとめて生成するの。必然的に for も使わなくて済むから。
> library(ggplot2)
> library(compiler)
>
> age <- data.frame(age = c(23, 25, 27, 27, 29, 26, 31, 31, 28, 29, 23, 26, 29,
+ 30, 30, 27, 26))
>
> getMean <- cmpfun(function(times, s = 10) {
+ sample_mean <- numeric(times)
+ for (i in 1:times) {
+ res <- sample(age$age, size = s, replace = TRUE)
+ sample_mean[i] <- mean(res)
+ }
+ return(sample_mean)
+ })
だからね,13秒もかかっちゃうわけよ。
> system.time(res1 <- getMean(1e+06))
ユーザ システム 経過
13.384 0.143 13.650
あいかわらず,ggplot の図も「きたない」し。
> ggplot(data.frame(res1), aes(x = res1)) + geom_histogram(binwidth = 0.1, colour = "white") +
+ geom_vline(x = mean(res1), color = "red") + annotate(x = 26, y = 40000,
+ geom = "text", label = paste(sep = ":", "mean", round(mean(res1), 1))) +
+ annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+ round(sd(res1), 3))) + xlim(24, 32)
まとめて標本を生成するとはこういうこと。
> getMean2 <- cmpfun(function(times, s = 10) {
+ res <- matrix(sample(age$age, size = s*times, replace = TRUE), times) # はい,ここですよ!
+ return(rowMeans(res)) # シミュレーション結果を返すのも一発で
+ })
で,こうなるわけさ。
> system.time(res2 <- getMean2(1e+06))
ユーザ システム 経過
0.231 0.057 0.287
> ggplot(data.frame(res2), aes(x = res2)) + geom_histogram(binwidth = 0.1, colour = "white") +
+ geom_vline(x = mean(res2), color = "red") + annotate(x = 26, y = 40000,
+ geom = "text", label = paste(sep = ":", "mean", round(mean(res2), 1))) +
+ annotate(x = 26, y = 30000, geom = "text", label = paste(sep = ":", "sd",
+ round(sd(res2), 3))) + xlim(24, 32)
定石なんだから,ああだこうだいわずに,いつでも採用すべき方策。
そもそも,この種の問題を解くのに,ブルート・フォースを使うのは良くない。
問題をよく考えれば,理論的な解がちゃんと存在するし,そこまで行かなくても,コンピュータ・パワーを利用すると瞬殺の答えが得られる。
以下のようなプログラムを書くと,n=10 なら 0.028 秒!!,n=398 でも 0.065 秒!!!
n=399 以上は,numeric が扱える上限を超える。gmp ライブラリあたりを使えば良い。
dice.n <- function(n=10) {
d <- c(1:6, 5:1)
for (n in 3:n) {
ld <- length(d)
d2 <- numeric(ld+5)
for (i in 1:6) {
d2[i:(i+ld-1)] <- d2[i:(i+ld-1)]+d
}
d <- d2
}
return(d)
}
> system.time({
+ z <- dice.n(10)
+ barplot(z)
+ })
ユーザ システム 経過
0.028 0.002 0.031
> system.time({
+ z <- dice.n(398)
+ barplot(z)
+ })
ユーザ システム 経過
0.065 0.004 0.068
もう正規分布だね。
「dice」にて,
> サイコロ3個ふって出目をカウント
> 6進法を使ってみた。
最後に
> サイコロ10個振ると?
とあるが,質問しているのか,挑戦しているのか,
少なくとも,掲載されているプログラムでは時間が掛かってやる気が起きない。
n <- 8
dice <- function(dec = 2, size = 2) (sum(mapply(as.integer, strsplit(dec2base(dec,
6, size), "")) + 1))
a <- rep(0, 6*n)
names(a) <- 1:(6*n)
for (i in 0:(6^n-1)) {
b <- dice(i, size = n)
a[b] <- a[b] + 1
}
print(a[n:(6*n)])
8個のサイコロを振ったら,126秒ほどかかる。10個振ると4500秒くらいもかかると思う。
=======
別のアプローチとして,expand.grid というのを使う。
x expand.grid(x, x, x) としてできるものは,各行は3つのサイコロの出目である。行数は 6^3 行である。
10個のサイコロを振るなら,expand.grid(x,x,x,x,x,x,x,x,x,x) とすれば,6^10行,10列の行列ができる。rowSums で出目の和をとり,table 関数で集計し,barplot で図を描く。
実行時間は 50 秒ほど。
> system.time({
+ n <- 10
+ x <- 1:6
+ y <- eval(parse(text=paste("expand.grid(", paste(rep("x", n), collapse=","), ")", collapse="")))
+ barplot(table(rowSums(y)), las=1)
+ })
ユーザ システム 経過
50.398 7.074 57.070
「mandelbrot set」に,以下のようなプログラムが掲載されている。
mandelbrot <- function(step = 0.1, pch = 19, col = 30, n = 2) {
a <- outer(seq(-2, 1, step), seq(-1, 1, step), function(x, y) complex(r = x,
i = y))
a <<- a[Mod(a) <= 2] #絶対値2以下のものだけにする(あまり意味はないが)
d <- c()
dc <- c()
for (c in a) {
z <- 0
for (i in 1:100) {
z <- z^n + c
if (Mod(z) >= 2) {
d <- c(d, c)
dc <- c(dc, i) #何回目で発散したか
break
}
}
}
# 回数で色分け
plot(d, pch = pch, xlab = paste("step:", step, "n = ", n), xlim = c(-2,
1), ylim = c(-1, 1), col = colors()[dc + col])
dc <<- dc
}
c() とは何事かと思うだろうが,NULL と同じである。名は体を表すのだから,NULL って書こうよ。
> identical(NULL, c())
[1] TRUE
それに,d <- NULL にしておき,後で d <- c(d, foo) などとドンドン要素を追加していくのは,プログラム作法に反している。計算時間が指数関数的に増えていく。
必要なだけのメモリを確保して,要素を代入していくのが定石だ。未だにこの辺りのことを体得していないプログラマがいるのは嘆かわしい。
それに,a, b, c, d, dc のように,意味のない変数名を使うのも困ったことだ。
<<- を使っているのも問題。必要ならば,list にして return 文で返すべきだ。
同でもよいことだけど,注釈記号の後には少なくとも空白1個置きましょう。気持ち悪いから。
ということで,一番目の指摘点だけを書き換えるだけで,プログラムの速度はかなり速くなる。
step により指数関数的に増えるのだが,step=0.01 のときは,元のプログラムの3倍ほどの速さになる。
mandelbrot <- function(step = 0.1, pch = 0, col = 30, n = 2) {
a <- outer(seq(-2, 1, step), seq(-1, 1, step), function(x, y) complex(real = x,
imaginary = y))
a <<- a[Mod(a) <= 2] # 絶対値2以下のものだけにする(あまり意味はないが)
len <- length(a) # 追加
d <- complex(len) # 書き直し
dc <- integer(len) # 書き直し
m <- 0 # 追加
for (c in a) {
z <- 0
for (i in 1:100) {
z <- z^n + c
if (Mod(z) >= 2) {
m <- m+1 # 追加
d[m] <- c # 書き直し
dc[m] <- i # 何回目で発散したか # 書き直し
break
}
}
}
# 回数で色分け
d <- d[1:m] # 書き直し
plot(d, pch = pch, xlab = paste("step:", step, "n = ", n), xlim = c(-2,
1), ylim = c(-1, 1), col = colors()[dc + col], asp=1, cex=10*step)
dc <<- dc[1:m]
}
「円を描く」にて
t <- seq(0, 2, pi/1000) * pi
x <- cos(t)
y <- sin(t)
plot(x, y, type = "l", asp = 1)
などと書いているが,
なんで,t <- seq(0, 2, pi/1000) * pi なのだろうか?意味不明である。
引数か何か誤解があるのでは?
それに,描かれた円は,閉じていない!!!
座標(1,0) のあたりを拡大すると,以下のようになっている!!
こんなことがないように,
t <- seq(0, 2*pi, length=1000)
でよいだろ
t <- seq(0, 2*pi, by=0.001)
は端点の処理という点で(t に 2*pi が含まれない),元のプログラムと同じで,今回の目的には合わない。
> t <- seq(0, 2, pi/1000) * pi
> sin(t[length(t)])
[1] -0.00611687 # 0に近くないとマズイ
> t <- seq(0, 2*pi, length=1000)
> sin(t[length(t)])
[1] -2.449294e-16 # ちゃんとやればこうなる。こうならないとだめ。