コイントスのシミュレーション

10 円硬貨を n 回投げたら，表は何回出るか？

断るまでもなく，使用する 10 円硬貨は裏表の出方に偏りはない，投げ方にも何らの癖もないという理想状態。

関数を書く

# n 回投げて，表の出た数を返す関数
cointoss(n) = count(isequal("head"), rand(["head", "tail"], n))   

cointoss (generic function with 2 methods)

# m 回の試行の結果（長さ m の整数ベクトル）を返す関数
simulation(m, n) = [cointoss(n) for i = 1:m]

simulation (generic function with 2 methods)

result = simulation(10000, 10);

using FreqTables
freq = freqtable(result)

11-element Named Vector{Int64}
Dim1  │ 
──────┼─────
0     │    9
1     │   98
2     │  423
3     │ 1187
4     │ 2074
5     │ 2407
6     │ 2057
7     │ 1194
8     │  434
9     │  106
10    │   11

# 結果を棒グラフで表示
using Plots
minx, maxx = extrema(result)
bar(names(freq), freq, xlabel="number of heads",　xticks=minx:maxx,
    ylabel="frequency", legend=false)
savefig("fig10.png")

一連の操作を若干の調整をして関数にまとめる。

using FreqTables, Plots
cointoss(n) = count(isequal("head"), rand(["head", "tail"], n))
simulation(m, n) = [cointoss(n) for i = 1:m]
function f(m, n)
    result = simulation(m, n);
    freq = freqtable(result)
    println(freq)
    minx, maxx = extrema(result)
    bar(names(freq), freq, xlabel="number of heads",
        ylabel="frequency", legend=false)
    savefig("fig$n.png")
end

f (generic function with 2 methods)

実行結果

f(10000, 2)

3-element Named Vector{Int64}
Dim1  │ 
──────┼─────
0     │ 2446
1     │ 5041
2     │ 2513

f(10000, 3)

4-element Named Vector{Int64}
Dim1  │ 
──────┼─────
0     │ 1265
1     │ 3740
2     │ 3702
3     │ 1293

f(10000, 5)

6-element Named Vector{Int64}
Dim1  │ 
──────┼─────
0     │  298
1     │ 1552
2     │ 3115
3     │ 3124
4     │ 1563
5     │  348

f(10000, 100)

38-element Named Vector{Int64}
Dim1  │ 
──────┼────
31    │   1
33    │   5
34    │   3
35    │   9
36    │  11
37    │  26
38    │  52
39    │  72
40    │ 111
41    │ 164
42    │ 231
43    │ 312
⋮         ⋮
58    │ 192
59    │ 157
60    │ 118
61    │  66
62    │  49
63    │  29
64    │  15
65    │  11
66    │   4
67    │   6
68    │   2
70    │   1

f(100000, 1000)

125-element Named Vector{Int64}
Dim1  │ 
──────┼───
429   │  1
439   │  2
440   │  1
441   │  2
442   │  1
443   │  3
444   │  7
445   │  6
446   │  8
447   │ 11
448   │  7
449   │ 14
⋮        ⋮
551   │ 14
552   │  9
553   │ 11
554   │  4
555   │  7
556   │  6
557   │  9
558   │  5
560   │  2
561   │  2
562   │  3
563   │  1

表の出る回数は二項分布にしたがう。コイントスの回数が多くなると，正規分布に近づいていく。

表の出る確率が 0.5 より偏っていても，コイントスの回数が多くなると，正規分布に近づいていく。

# 表の出る確率が 0 < p < 1 のコインを n 回投げて，表の出た数を返す関数
# 前述の関数と同名であるが，引数が 2 個の場合。実装法を変えたが同じ機能を持つ。
cointoss(n, p) = sum(rand(n) .< p)

cointoss (generic function with 2 methods)

# m 回の試行の結果（長さ m の整数ベクトル）を返す関数
simulation(m, n, p) = [cointoss(n, p) for i = 1:m]

simulation (generic function with 2 methods)

function f(m, n, p=0.5)
    result = simulation(m, n, p);
    freq = freqtable(result)
    println(freq)
    minx, maxx = extrema(result)
    bar(names(freq), freq, xlabel="number of heads",
        ylabel="frequency", legend=false)
    savefig("fig2-$n.png")
end

f (generic function with 2 methods)

f(1000, 10, 0.2)

8-element Named Vector{Int64}
Dim1  │ 
──────┼────
0     │  83
1     │ 288
2     │ 291
3     │ 222
4     │  85
5     │  23
6     │   7
7     │   1

f(10000, 100, 0.2)

31-element Named Vector{Int64}
Dim1  │ 
──────┼────
6     │   2
7     │   1
8     │   4
9     │  14
10    │  35
11    │  75
12    │ 107
13    │ 222
14    │ 340
15    │ 526
16    │ 621
17    │ 804
⋮         ⋮
25    │ 443
26    │ 319
27    │ 213
28    │ 141
29    │  85
30    │  49
31    │  27
32    │  24
33    │   6
34    │   7
35    │   1
36    │   4