『確率4万分の1、県の入札くじに6回連続当選「奇跡的」』なのかどうかを調べよう
動機
面白いニュースを見かけたので計算してみましたなポストです。記事のこの部分に注目しました。
6回のくじにはそれぞれ制限価格で入札した3~8の鑑定業者が参加した。6回連続で当たる確率を計算すると約4万分の1となる。
奇跡ってどのくらいの割合で起きるのかな?もうすぐクリスマスだしってことでRでやってみます。
まず引用箇所を表にすると、
i回目 | 鑑定業者数 X_i |
---|---|
1 | X_1 |
2 | X_2 |
3 | X_3 |
4 | X_4 |
5 | X_5 |
6 | X_6 |
ここでX_i ∈ [3,8]の整数で、「6回連続で当たる確率を計算すると約4万分の1」を満たすX_iであって欲しいってことでX_iを全部掛け算したら4万になれば良いでしょう。
とりあえずそのような鑑定業者数の組み合わせを力技で探すと、
kgyousya <- c(3,4,5,6,7,8) kaisuu <- data.frame( x1=kgyousya ,x2=kgyousya ,x3=kgyousya ,x4=kgyousya ,x5=kgyousya ,x6=kgyousya ) kaisuu.d <- expand.grid(kaisuu) kaisuu.m <- as.matrix(kaisuu.d) kaisuu.prod <- apply(kaisuu.m,MARGIN = 1,prod) table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] ) #> table( kaisuu.prod[(kaisuu.prod > 39000) & (kaisuu.prod < 41000)] ) # #39200 40000 40320 40960 #180 15 720 60
915ケースありました。
次に、掛け算して39,200〜40,960となるような鑑定業者数の組み合わせ:915ケースで、 各ケースを10000回シミュレーションする。 そして各ケースで確率4万分の1以下が10000回中何回起きるかを調べる。(←奇跡の割合) 最後に奇跡の割合の分布を観察しましょう。
# 候補のケース数: 915 #> nrow(kaisuu.sub) #[1] 915 mycase <- nrow(kaisuu.sub) results <- numeric(mycase) # シミュレーション回数 ntimes <- 1:10000 result <- numeric( length(ntimes) ) for(j in 1:mycase){ for(i in ntimes){ np <- sample(c(0,1),size=6,replace = TRUE, prob=c(.5,.5)) result[i] <- prod( abs( np - (1/kaisuu.sub[j,]) ) ) } results[j] <- sum( result < 1/40000 )/length(result) }
結果
度数分布表
> table(results) results 0 0.005 0.006 0.007 0.008 0.009 0.01 0.011 0.012 0.013 0.014 0.015 0.016 0.017 0.018 180 1 3 7 8 17 21 56 63 66 64 54 74 74 49 0.019 0.02 0.021 0.022 0.023 0.024 0.025 0.026 0.027 0.028 0.029 46 45 30 18 19 7 4 4 3 1 1
要約値
> summary(results) Min. 1st Qu. Median Mean 3rd Qu. Max. 0.00000 0.01000 0.01400 0.01223 0.01700 0.02900
どんなケースでも平均すれば1%ぐらいは起きるのでは?という見方で良いのかなぁ・・・。
プロット
感想
こういうの式だけ立てて計算できるようになりたいですなぁ。
2018年11月 統計検定2級 問9
この問題をなぜ解くのか
前にこんなこと
統計検定であって確率検定ではないので、サイコロを7回投げる時2以下の目が出る確率を解けるようになってもなんだかなぁ
統計検定2級に合格したので勉強法やら参考書などを書いておきます。 - axjack's blog
を書きました。がしかし、統計検定 準1級の勉強をしていく中「やっぱ確率もちゃんと解けないとダメじゃん😢」と気持ちが一転したので改めてちゃんと解いてみようと思った次第です。
問題文の概要
どの面も等確率で出る6面のサイコロを7回投げる。2以下の目が出る回数をXとする。
- 問(1):P(X = x+1)とP(X = x)の比、i.e. を求めよ
- 問(2):P(X = x) が最大となるxは?
解答
問(1)
を出目が2以下の時に1、そうでない時を0とするとはベルヌーイ分布:Bin(1,p)に従う。ここでp=1/3である。この時、 は二項分布:Bin(n,p)に従う。よって、確率質量関数はとなる。ここで、n=7, p=1/3を代入するとP(X=x)とP(X=x+1)はそれぞれ、
となる。すると、
となる。よって、
問(2)
P(X=x)が最大となるxを を用いて考えてみると、Q(x)が1より大であればP(X=x+1) > P(x)であると言えるし1未満であればP(X=x+1) < P(x)と言える。具体的にxへ0から6を代入すると、
x | Q(x) | 判定 |
0 | 7/2 > 1 | P(1) > P(0) |
1 | 6/4 > 1 | P(2) > P(1) |
2 | 5/6 < 1 | P(3) < P(2) |
3 | 4/8 < 1 | P(4) < P(3) |
4 | 3/10 < 1 | P(5) < P(4) |
5 | 2/12 < 1 | P(6) < P(5) |
6 | 1/14 < 1 | P(7) < P(6) |
となる。関数の増減表のように判定の列を眺めれば、P(x)はP(2)がピークであることからP(X = x) が最大となるxは2である。
超幾何分布
何度やっても忘れるのでブログに書いて覚えよう。
超幾何分布の確率質量関数
計算例
その1
男50人女50人から10人を選ぶ。10人のうち男3人女7人となる確率pは、超幾何分布を用いて
その2
統計検定準一級2017年問10より引用。表の80人から30人を無作為抽出する。その30人のうち男性かつ就職している人数Xは超幾何分布に従う。
\ | 就職 | 非就職 | 計 |
男性 | 38 | 3 | 41 |
女性 | 30 | 9 | 39 |
計 | 68 | 12 | 80 |
- 80人から30人を選ぶ
- 「男性かつ就職」である人:38
- 「男性かつ就職」でない人:80-38 = 42
と整理して、
となる。
その3
同じ表を用いて、80人から20人を無作為抽出する。その20人のうち「女性かつ非就職」な人数Xは超幾何分布に従う。
- 80人から20人を選ぶ
- 「女性かつ非就職」である人:9
- 「女性かつ非就職」でない人:80-9 = 71
と整理して、
となる。
Rで計算する
choose関数で立てた式が超幾何分布の質量関数dhyperと同じことを確認。
# その1 > choose(50,3)*choose(50,7)/choose(100,10) [1] 0.1130964 > dhyper(3,50,50,10) [1] 0.1130964
dhpyer(x,m,n,k)の引数*1は、
- x: vector of quantiles representing the number of white balls drawn without replacement from an urn which contains both black and white balls.
- m: the number of white balls in the urn.
- n: the number of black balls in the urn.
- k:the number of balls drawn from the urn.
連の検定(Runs Test for Detecting Non-randomness)
はじめに
連の検定についてよく分からなかったのでRで実装して確かめてみました。
連とは?
2値{A,B}を取る系列ABAABBAAABBBAB
があった時、これをA | B | AA | BB | AAA | BBB | A | B
と連続する同じ値ごとに分割できるように見える。この時同じ文字または記号のひと続きを連(run)
*1と言い、またひと続きの個数を連の長さ
*2と言う。この場合、Aの連が4つ(= {A,AA,AAA,A})とBの連が4つ(= {B,BB,BBB,B})あることが分かる。AAAAAAAAABBBBBBBBBBB
であれば、Aの連が1つでBの連が1つとなる。
確かめたこと
Aの数(= Na)とBの数(= Nb)が与えられた時、N = Na+Nbの大きさを持つ系列は、
を確かめた。
ソース
# 連の検定#### rm(list=ls()) # -------------------- # 関数: generateRun #### # -------------------- # 引数 # Na_: Aの個数 # Nb_: Bの個数 # 戻り値 # 連の数を返す関数 # -------------------- generateRun <- function(Na_, Nb_){ return( function(dummy){ N <- Na_ + Nb_ Arep <- rep(TRUE,Na_) Brep <- rep(FALSE,Nb_) AB <- c(Arep,Brep) xx <- sample(AB, N, replace = FALSE) yy <- xx[1:(N-1)] zz <- xx[2:N] return( sum( yy != zz) + 1 ) } ) } # シミュレーション #### ## 各種変数設定 Na <- 9 # Aの個数 Nb <- 11 # Bの個数 sim_N <- 10000 # シミュレーション回数 ## 平均値と分散 ( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値 ( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散 # シミュレーション実行 Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) ) # 要約 table(Result.Run) summary( Result.Run ) var( Result.Run ) # グラフ表示 hist(Result.Run,breaks = seq(min(Result.Run),max(Result.Run),1)) par(new=T) # 正規分布を重ね合わせ curve(dnorm(x,t_mean,sqrt(t_var)),min(Result.Run),max(Result.Run),xlab="",ylab="",xaxt="n", yaxt="n",col="red",bty="n")
結果の確認
数値の確認
> # シミュレーション #### > ## 各種変数設定 > Na <- 9 # Aの個数 > Nb <- 11 # Bの個数 > sim_N <- 10000 # シミュレーション回数 > ## 平均値と分散 > ( t_mean <- (2*Na*Nb)/(Na+Nb) + 1 ) # 平均値 [1] 10.9 > ( t_var <- (2*Na*Nb*(2*Na*Nb - Na - Nb))/((Na+Nb)^2 * (Na+Nb-1)) ) # 分散 [1] 4.637368 > # シミュレーション実行 > Result.Run <- sapply( 1:sim_N, generateRun(Na,Nb) ) > # 要約 > table(Result.Run) Result.Run 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 1 3 13 28 167 312 852 1218 1725 1699 1699 1144 704 279 121 30 5 > summary( Result.Run ) Min. 1st Qu. Median Mean 3rd Qu. Max. 2.00 9.00 11.00 10.89 12.00 18.00 > var( Result.Run ) [1] 4.666794
プロット
結果の考察
- 計算した平均値・分散がシミュレーションの結果とだいたい合っていることがわかった。
- ヒストグラムとカーブプロットが大まかにあっている。
おわりに
帰無仮説の下での平均値・分散の式が奇妙で覚えられないし式を見てもあまり実感が湧かなかったものの、シミュレーションしてみると確かに平均値・分散がそれっぽく合うのでこの式の妥当性を確認することができた。疑わしい系列があってもこれで確かめることができそうだ。
参考
連の検定関連
R言語のplot関連
- 49. 図の重ね合わせ
- 図の重ね合わせ方法を参考にしました。
- 正規分布のグラフをRで描く。【curve()の使い方】
- 正規分布をcurveで書く方法を参考にしました。
- R 使い方 軸・ラベルの調整(向き・サイズ・色など) グラフの描き方
- plotを重ねる時に軸やラベルを消す方法を参考にしました。
自由度調整済み決定係数を決定係数で表す
やること
タイトルの通り式変形をするだけです。
式変形
2019年8月をワールドワイドに考える
はじめに
UTCを考慮すると途端にわけわからなくなったので整理する。
UTCとは
wikiで協定世界時 - Wikipediaをみませう。
UTCってプラマイいくつまであるの?
時間帯 (標準時) - WikipediaによるとからUTC-12からUTC+14まである。日本はUTC+9
2019年8月をUTCで列挙する
以下の時間帯で考える。
すると次の表のようになる
UTC-12 | UTC0 | UTC+9 | UTC+14 |
---|---|---|---|
7/30 22:00 | 7/31 10:00 | 7/31 19:00 | 8/01 00:00 |
7/31 03:00 | 7/31 15:00 | 8/01 00:00 | 8/01 05:00 |
7/31 03:00 | 8/01 00:00 | 8/01 09:00 | 8/01 14:00 |
8/01 00:00 | 8/01 12:00 | 8/01 21:00 | 8/02 02:00 |
… | … | … | … |
8/30 22:00 | 8/31 10:00 | 8/31 19:00 | 9/01 00:00 |
8/31 03:00 | 8/31 15:00 | 9/01 00:00 | 9/01 05:00 |
8/31 03:00 | 9/01 00:00 | 9/01 09:00 | 9/01 14:00 |
9/01 00:00 | 9/01 12:00 | 9/01 21:00 | 9/02 02:00 |
結局2019年8月とは
- UTC+14が世界で最も早く8/1を迎えるので、UTC0で考えれば7/31 10:00(⇔日本時間 7/31 19:00)
- UTC-12が世界で最も遅く9/1を迎えるので、UTC0で考えれば9/1 12:00(⇔日本時間 9/1 21:00)
の時間帯を抽出すればとりあえず取りこぼすことはなさそうだ。