Livebook で 確率分布を描く (Elixir, VegaLite)

Elixir に **Statistics** という統計のライブラリがあるので試してみました。２項分布と正規分布という有名な確率分布を計算して、**Livebook** 上の **VegaLite** で描画してみる試みです。
[Statistics](https://hexdocs.pm/statistics/Statistics.html)

確率変数については以下の過去記事を参照してください。
[Juliaで学ぶ確率変数(1) - 確率変数の定義](https://qiita.com/sand/items/afb5976ad28120d216c4)

当然、この分野では **Python** や **Julia** の充実には、まだまだ届いていない感じです。
[Elixir and Machine Learning: Nx v0.1 released!](https://dashbit.co/blog/elixir-and-machine-learning-nx-v0.1)



**関連記事**
[VegaLite で埼玉県を切る( Elixir, Livebook) - Qiita](https://qiita.com/sand/items/a8b99ee8b68f49b5784a)
[Vega-Lite の View Composition (Elixir, Livebook) - Qiita](https://qiita.com/sand/items/f0e11ed907287be6436a)
[VegaLite の基礎 (Elixir, Livebook) -Qiita](https://qiita.com/sand/items/3433061b1966d6f8edb2)
[Livebook 事始め (Elixir) - Qiita](https://qiita.com/sand/items/6753c1782774b45daf6b)


まずは必要なインストールを行っておきます。
```elixir:
Mix.install([
  {:vega_lite, "~> 0.1.6"},
  {:kino_vega_lite, "~> 0.1.4"},
  {:jason, "~> 1.2"},
  { :statistics, "~> 0.6.2"}
])

alias VegaLite, as: Vl
```


# 1. ２項分布　Statistics.Distributions.Binomial

[Juliaで学ぶ確率変数(2) - ２項分布（離散型）](https://qiita.com/sand/items/9eda536b4e8f9bad0055)

**cdf(n, p)** - The cumulative density function （確率密度関数）
**pmf(n, p)** - The probability mass function （確率関数）

## 1-1. ２項分布 例題１
::: note info
左利きの人が1%いる集団から1人選び、右利きか左利きかをみる。それを200回繰り返したときに4人以上が左利きである確率を計算する。
:::

* n = 200
* p = 0.01

とりあえず分布図を描いてみます。

```elixir:
n = 200
p = 0.01
xs = 1..20
ys =
  xs
  |> Enum.map(&Statistics.Distributions.Binomial.pmf(n, p).(&1))

Vl.new()
|> Vl.data_from_values(x: xs, y: ys)
|> Vl.mark(:line)
|> Vl.encode_field(:x, "x", type: :quantitative )
|> Vl.encode_field(:y, "y", type: :quantitative )
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61915/4c28a8a3-6c38-635e-6355-2af271cc7797.png)

4人の時の値だけを計算させます。

```elixir:
Statistics.Distributions.Binomial.pmf(n, p).(4)
```

```
0.09021970192446989
```


## 1-2. ２項分布 例題２

::: note info
ある工場で製造される工業製品には、0.1%の割合で不良品が混ざる。この工業製品1万個の中の不良品の個数をXであらわす。1万個の製品は独立に製造されたと考えるとき、Xの平均と分散を求めよ。
:::

* n = 10000
* p = 0.001

分布図を描きます。

```elixir:
n = 10000
p = 0.001
xs = 1..30
ys =
  xs
  |> Enum.map(&Statistics.Distributions.Binomial.pmf(n, p).(&1))

Vl.new()
|> Vl.data_from_values(x: xs, y: ys)
|> Vl.mark(:line)
|> Vl.encode_field(:x, "x", type: :quantitative )
|> Vl.encode_field(:y, "y", type: :quantitative )
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61915/77590c4e-a8ad-15e3-1a4d-f81301a6c50e.png)


平均値 = n * p、と 分散 = n * p * (1-p) を求めます。
```elixir:
mean = n*p
var = n*p*(1-p)
[mean, var]
```

```
[10.0, 9.99]
```

# 2. 正規分布 Statistics.Distributions.Normal

[Juliaで学ぶ確率変数(5) - 正規分布（連続型）](https://qiita.com/sand/items/4c3ef0b0504704074681)

* **cdf(mu, sigma)** - The cumulative density function（確率密度関数）
* **pdf(mu, sigma)** - Probability density function（確率密度関数）


::: note info
日本の成人男性の平均身長 μ＝171cm、標準偏差 σ＝6cmとする。
正規分布に従うと仮定した場合、日本中からランダムに1人選ばれた成人男性の身長が165cm以上171cm以下である確率は何%か？
:::



* μ = 171
* σ ＝ 6


分布図を描きます。

```elixir:
mu = 171
sigma = 6
xs = 150..200
ys =
  xs
  |> Enum.map(&Statistics.Distributions.Normal.pdf(mu, sigma).(&1))

Vl.new()
|> Vl.data_from_values(x: xs, y: ys)
|> Vl.mark(:line)
|> Vl.encode_field(:x, "x", type: :quantitative )
|> Vl.encode_field(:y, "y", type: :quantitative )
```

![image.png](https://qiita-image-store.s3.ap-northeast-1.amazonaws.com/0/61915/bef62f36-783b-6331-8d99-d3b25c1fd355.png)


165cm以上171cm以下である確率は、x=171 と x=165 の間の面積となるから、単に引き算を行えばよい。

```elixir:
mu = 171
sigma = 6
Statistics.Distributions.Normal.cdf(mu, sigma).(171) -
  Statistics.Distributions.Normal.cdf(mu, sigma).(165) 
```

```
0.3413447366676363
```


今回は以上です。