トップ 履歴 一覧 ソース 検索 ヘルプ ログイン

R 3.6以前と以後の乱数生成器の動作を比べてみる


キーワード

最終更新時間:2021年05月14日 01時13分45秒
アフィリエイト・広告について
プライバシーポリシー

はじめに

1、2年前に、Rのバージョンアップに合わせて、「Rの乱数生成のアルゴリズムが変わった」と、少し話題になりました。その時は、「へー、まぁよりランダムに、いい感じになったんでしょ」というくらいにしか捉えていませんでしたが、あらためて、どういうことか調べてみました。といっても、私は数学はまるでできないので、乱数生成の理論については触れません。ただ、自分で調べたときに、検証した結果のグラフなどが出てこなかったので、興味がある方の参考になれば、ということで。

上記の文章を、もう少し正確に記述すると、2019年4月にリリースされた、R 3.6.0において、離散一様分布を生成するアルゴリズムが変更されました。リリースのアナウンスには、次のようにあります。

CHANGES IN R 3.6.0:

SIGNIFICANT USER-VISIBLE CHANGES:

* The default method for generating from a discrete uniform

distribution (used in sample(), for instance) has been changed.

This addresses the fact, pointed out by Ottoboni and Stark, that

the previous method made sample() noticeably non-uniform on large

populations. See PR#17494 for a discussion. The previous method

can be requested using RNGkind() or RNGversion() if necessary for

reproduction of old results. Thanks to Duncan Murdoch for

contributing the patch and Gabe Becker for further assistance.

The output of RNGkind() has been changed to also return the

'kind' used by sample().

https://stat.ethz.ch/pipermail/r-announce/2019/000641.html

ざっと要約すると、sample() 関数などで使われている、これまで (3.6以前) の乱数生成アルゴリズムは、母集団が巨大な場合において、本当に一様ではなかった、とあります。また、以前のアルゴリズムも、指定すれば使用できる、とあります。そこで、変わってからだいぶ経ちますが、3.6以前と以降で、乱数生成の結果がどう変わるのかを確認してみます。

Rのバージョンの確認

まず、現在のRのバージョンを確認しておきます。


私の実行環境では、以下のようになります。

platform       x86_64-pc-linux-gnu         
arch           x86_64                      
os             linux-gnu                   
system         x86_64, linux-gnu           
status                                     
major          4                           
minor          0.5                         
year           2021                        
month          03                          
day            31                          
svn rev        80133                       
language       R                           
version.string R version 4.0.5 (2021-03-31)
nickname       Shake and Throw           

乱数生成器の確認

このバージョン (4.0.5) における、乱数生成器の設定を確認します。

実行結果は、以下のようになります。

[1] "Mersenne-Twister"   "Inversion"   "Rounding"

RNGkind() 関数は、現在設定されている、(1) 連続一様乱数を生成するアルゴリズム、(2) 正規乱数を生成するアルゴリズム、(3) 離散一様乱数を生成するアルゴリズムを返します。つまり、R 4.0.5では、(1) 連続一様乱数の生成にMersenne-Twister法、(2) 正規乱数にInversion法、そして (3) 離散一様乱数にRejection法を使っています。(1) と (2) については、昔から特に変更はないので、ここでは (3) のRejection法による乱数と、それ以前の方法 (Rounding) による乱数を比較してみます。

現在のRにおける乱数生成と結果

ここでは、乱数生成の方法を変更した、というパッチに記載されている方法で検証してみます。

コードでは、5で割り切れる巨大な数 m を生成しています。ここでは、1342177280になりました。そして、sample() 関数で、1から m までの数値から100万個をサンプリングします。replace = TRUE なので、同じ値が複数回サンプリングされる可能性があります。次に、その乱数列を5で割った余りを求めて、table() 関数で集計します。理論上、0から4までが均等に得られるはずです。結果は以下のようになります。(乱数なので、お手元で実行すると結果は変わります)

     0      1      2      3      4 
200029 200038 200458 200245 199230 

R 4.0.5のデフォルトの乱数生成アルゴリズム (Rejection) で乱数を生成した結果

最も多い値と、最も少ない値の差は0.6%程度と、確かに一様であることがわかります。

R 3.6以前における乱数生成と結果

では、R 3.6以前の方法 (Rounding) ではどうなるでしょうか。上述のように、以前の方法も、指定すれば使用できます。

実行すると、以下のように警告が出ます。

Warning message:
In RNGkind(sample.kind = "Rounding") : non-uniform 'Rounding' sampler used

この状態で、同じく sample() 関数で、1から m までの数値から100万個をサンプリングし、得られた乱数を5で割って、結果を集計します。

結果は以下のようになります。

     0      1      2      3      4 
187262 250030 187623 187840 187245

R 3.6以前の乱数生成アルゴリズム (Rounding) で乱数を生成した結果

明らかに、一様ではありません。最も多い値と、最も少ない値の差が、33%あまりに上っています。

なぜそうなるのか、という理論的なところは私にはわかりませんが、事実として確かに、R 3.6以前の乱数生成アルゴリズムには問題があることがわかりました。


カテゴリ: [R, データサイエンス, 統計]