numpy.randomのGeneratorをためしてみる

初投稿です。

NumPy v1.18が2020年5月末にリリースされました。それで改めてドキュメントを読んでいて気づいたのですが、いつのまにか、統計/機械学習でお世話になっている`numpy.random`モジュールに新たに`Generator`なるクラスが追加され、これまで馴染んでいた`numpy.random.random`などの関数が非推奨となったようです。

しかも、この変更はNumPy1.17で行われていたようです。今更気づく情弱......

https://numpy.org/doc/1.17/reference/random/index.html

個人的には結構大きな変更だと思ったので、メモがてら何が変わったのか要点だけでも押さえておきたいと思います。

以下の記事は、NumPy v1.18のドキュメントに基づきます。

## 新しいランダムサンプリング

NumPyのドキュメントによると、これまで

```python
from numpy.random import random

random(100)
```

と書いていたランダムサンプリングのコードのかわりに、v1.17以降では

```python
from numpy.random import default_rng

rg = default_rng()
rg.random(100)
```

と書くことが推奨されています。
なんだかオブジェクト指向っぽいですね。

ちなみにシードを設定する場合はこうなります。

```python
from numpy.random default_rng

rg = default_rng(123)
rg.random(100)
```

これに対応する古いコードはこうですね。

```python
from numpy.random import RandomState, random, seed

rg = RandomState(123)
rg.random(100)

# グローバルシードを設定する場合
seed(123)
random(100)
```

## 新しい方法のメリット

ぶっちゃけこれ見たときに「行数減らへんやんけ」「`default_rng`ってなんやねん」「おどれなんのためにこんなもん使わなあかんのや」と思いましたが、なかなかのメリットがあるみたいです。

### メリット1: 速い

まず、速くて効率的です。私のPCの環境だと、

```python
In [1]: from numpy.random import standard_normal, default_rng

In [2]: rg = default_rng()

In [3]: %%timeit
   ...: rg.standard_normal(1000000)  # newer 
   ...:  
   ...:                                                                         
24.8 ms ± 709 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

In [4]: %%timeit
   ...: standard_normal(1000000)  # older
   ...:                                                                         
52.4 ms ± 835 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)
```

という結果になりました。なんと実行時間が半分になってます!

この理由は、内部で使っている乱数生成器が変わったからだそうです。v1.16以前では乱数生成器にメルセンヌ・ツイスタを利用していましたが、v1.17以降ではPCGという乱数生成器を利用しています。

PCGは2014年に発表された乱数生成器で、速くて効率的で性質のよい疑似乱数の生成が可能なんだそうです。私は詳しくないので、気になる方は[Wikipedia](https://ja.wikipedia.org/wiki/Permuted_congruential_generator)や[オフィシャルページ](https://www.pcg-random.org/)などで調べてみてください。解説qiita記事待ってます。

### メリット2: 乱数生成器が選べる

v1.16以前では乱数生成器にメルセンヌ・ツイスタ以外は使えませんでした。そこには選択肢なんてありませんでした。もしかしたらC言語が使える魔術師ならなんとかできたのかもしれませんが、わたしのように愚かで哀れでC言語を完全に理解することのできない1人のエンドユーザーには、与えられたメルセンヌ・ツイスタという道具を利用することしかできなかったのです。

しかし、v1.17以降は違います。

v1.17以降、私たちには乱数生成器を選ぶ権利が与えられました。私たちは、デフォルトのPCGだけでなく、懐かしのメルセンヌ・ツイスタを使うこともできるし、他にもいくつかの乱数生成器を選ぶことができます。

そのキモとなる仕組みが、`Generator`と`BitGenerator`です。

## `(Bit)Generator` is 何

`BitGenerator`クラスは、PCGやメルセンヌ・ツイスタのような乱数生成器の抽象クラスのようなもの、みたいです(フワフワ理解)。NumPyでは、`PCG64`、`MT19937`、`Philox`、`SFC64`といった疑似乱数生成器が`BitGenerator`のサブクラスのように使えます。

試してはいませんが、おそらく`BitGenerator`クラスを継承することで、カスタム乱数生成器を作成して使うこともできるはずです。やったねPythonちゃん!乱数が増えるよ!

そして`Generator`クラスは、この`BitGenerator`を一様分布や正規分布などの特定の確率分布に結びつける役割を担っています。

`random`や`normal`などの馴染み深い関数は`Generator`クラスのメソッドに生えています。`gumbel`や`weibull`のようなちょっとマニアックな確率分布の関数も生えています。嬉しいですね。

では、実際に`Generator`を使って乱数を作ってみましょう。

```python
from numpy.random import Generator, PCG64, MT19937

# 乱数生成器にPCGを使う
rg_pcg = Generator(PCG64())
rg_pcg.random(100)

# 乱数生成器にメルセンヌ・ツイスタを使う
rg_mt = Generator(MT19937())
rg_mt.random(100)
```

### `default_rng` is 何

最初の方に出てきた`default_rng`は、乱数生成器に`PCG64`を使う`Generator`のコンストラクタです。

IPythonでドキュメントを見てみましょう。

```python
In [1]: default_rng?
Docstring:
Construct a new Generator with the default BitGenerator (PCG64).
```

そういうわけで、以下のコードは実質同じということになります。

```python
# default_rng版(seed有り)
rg = default_rng(123)

# Generator版(seed有り)
rg = Generator(PCG64(123))
```

情弱ゆえv1.17がリリースされてから1年遅れのキャッチアップとなりましたが、少しでもお役に立てれば幸いです。

良いPythonライフを!