Edited at

Juliaで学ぶ統計学入門(第4章:確率)

More than 5 years have passed since last update.


はじめに

Facebookグループ『東大の統計学教科書を読み進める会』で、

統計学入門(通称:赤本)の読み進めを行っています。

赤本の処理をJuliaでどう実装するかが気になったので、コードを書き連ねています。


統計学入門(通称:赤本)とは?

「東京大学教養学部統計学教室 編」の統計学入門書です。古典です。

参考:統計学入門


Juliaとは?

科学計算向けのプログラミング言語です。ポストR・ベターMATLABという扱いです。

参考:Juliaの公式サイト


調査対象の章

第4章:確率


コード


sec04.jl

Pkg.add("Distributions")

using Distributions
Pkg.add("StatsBase")
using StatsBase

# 確率
ary=[1,1,1,1,2,2,2,3,3,4]
proportions(ary,minimum(ary):maximum(ary))
# =>
# 4-element Array{Float64,1}:
# 0.4
# 0.3
# 0.2
# 0.1

# 乱数
n = 3
rand() # 1個の乱数(一様分布)、 0 <= rand() < 1
rand(n) # n個の乱数(一様分布)
randn() # 1個の乱数(正規分布)、 0 <= randn() < 1
randn(n) # n個の乱数(正規分布)

# 切り上げ・切り捨て
ceil(rand(),n) # 小数第n位に切り上げ
floor(rand(),n) # 小数第n位に切り捨て

# 階乗(factorial)
factorial(n) # n!

# 組み合わせ(combination)数
n=5; r=2;
binomial(n, r) # 組み合わせ数、nCr

# 順列(permutation)数
n=5; r=2;
binomial(n, r) * factorial(r) # 順列数、nPr

# 円周率
pi # => 3.1415926535897...

# 練習問題

## 4.1 ド・メレの問題
### ⅰ)
d = Binomial(4, 1/6)
pdf(d, 0) # => 0.4822530864197532
sum(pdf(d, 1:4)) # => 0.5177469135802468

### ⅱ)
d = Binomial(24, 1/36)
pdf(d, 0) # => 0.5085961238690967
sum(pdf(d, 1:36)) # => 0.4914038761309033

## 4.2 ホイヘンスの14問題の第11
chs = repeat(collect([1:6]), inner=[6]) + repeat(collect([1:6]), outer=[6])
n = size(chs, 1)
cnt = count(x -> (x == 12), chs)
prp = cnt/n
i = 1
while !(0.9 <= (1 - pdf(Binomial(i, prp), 0)))
i += 1
end
print(i) # => 82

## 4.3
binomial(30, 15) # => 155117520

## 4.4 誕生日の問題
p(r) = 1 - reduce((prd,x) -> prd * (1 - (x/365)), 1, 1:(r-1))
# 1 - ((365-1)/365 * (365-2)/365 * ... * (365-(r-1)/365))
rs = vcat([5:5:15], [20:1:25], [30:5:40], [50,60])
map(p, rs)
# =>
# 14-element Array{Float64,1}:
# 0.0271356
# 0.116948
# 0.252901
# 0.411438
# 0.443688
# 0.475695
# 0.507297
# 0.538344
# 0.5687
# 0.706316
# 0.814383
# 0.891232
# 0.970374
# 0.994123

## 4.5 行動科学と確率
d = Binomial(45, 1/2) #向かいの席or隣の席の2択と解釈する場合
pdf(d, 30) # => 0.009801721761959955 (1%有意水準で有意)
d = Binomial(45, 2/3) #向かいの席or右隣の席or左隣の席の3択と解釈する場合
pdf(d, 30) # => 0.12534170054634222 (5%有意水準で有意でない)

## 4.6 ガリレイの問題
chs = repeat(collect([1:6]), inner=[6 * 6]) + repeat(collect([1:6]), outer=[6 * 6]) + repeat(repeat(collect([1:6]), inner=[6]), outer=[6])
cnt = count(x -> (x == 9), chs) # => 25 ([x,x,x]と[x,x,y]は組み合わせ数が違う。[x,y,x][y,x,x]のケースがある)
cnt = count(x -> (x == 10), chs) # => 27

## 4.7 ガンの診断
## ⅰ)
(0.005 * 0.95) / (0.005 * 0.95 + 0.995 * 0.05) # => 0.0871559633027523
#P(C|A) = P(A*C) / P(A)
#P(A) = P(A*C) + P(A*!C)
#P(A*C) = P(C) * P(A|C)
#P(A*!C) = P(!C) * P(A|!C)0.90 * 0.995 *
#P(AC) = 1 - P(!A|!C)
#P(C|A) = P(C) * P(A|C) / (P(C) * P(A|C) + P(!C) * (1 - P(!A|!C)))
## ⅱ)
- (0.90 * 0.995) / (0.90 * 0.005 - 0.90 * 0.995 - 0.005) # => 0.9994419642857142
#P(C|A) = P(C) * P(A|C) / (P(C) * P(A|C) + P(!C) * (1 - P(!A|!C)))
#0.90 = 0.005 * P(A|C) / (0.005 * P(A|C) + 0.995 * (1 - P(A|C)))
#P(A|C) = - (0.90 * 0.995) / (0.90 * 0.005 - 0.90 * 0.995 - 0.005)



感想

JuliaLangで調べた事をQiitaにまとめるよりも、

メモ代わりのJLJpWikiを作ってそこに記載する方が良い気がしてきた。