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

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

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

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

## 関数を書く



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




    cointoss (generic function with 2 methods)




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




    simulation (generic function with 2 methods)




```julia
result = simulation(10000, 10);
```


```julia
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




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

![fig10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/fe21768c-30bb-88ab-d360-fe49fe83f92c.png)

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


```julia
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)



## 実行結果


```julia
f(10000, 2)
```

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

![fig2.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/b4ee0e78-dabe-ac17-27b9-cf465546bca1.png)


```julia
f(10000, 3)
```

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

![fig3.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/b492d732-f753-2bb5-e974-3c09a78d5a6f.png)


```julia
f(10000, 5)
```

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

![fig5.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/5415bccc-d1e9-0f69-e662-75857b49a04d.png)


```julia
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

![fig100.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/d555917a-88dd-b5a0-4cf0-b1afc5053e38.png)


```julia
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

![fig1000.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/305496c9-c77a-609b-7099-8c85676b6d8e.png)



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

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


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




    cointoss (generic function with 2 methods)




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




    simulation (generic function with 2 methods)




```julia
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)




```julia
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

![fig2-10.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/086676f0-7153-cc7a-a396-aa94e4158c46.png)


```julia
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



![fig2-100.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/1653782/d3db52ae-713d-386f-b692-018cd183e22e.png)