Juliaで学ぶ統計学入門（第4章:確率）

はじめに

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)  # ｎ個の乱数（一様分布）
randn()  # 1個の乱数（正規分布）、 0 <= randn() < 1
randn(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を作ってそこに記載する方が良い気がしてきた。