洋食の日記

「だ・である」調ではなく「です・ます」調で書きはじめれば良かったなと後悔してる人のブログです

確率分布に従う乱数によるNumo::NArrayを生成するNumo::Randomを作った

はじめに

正規分布や二項分布などの確率分布に従う(疑似)乱数によるNumo::NArrayを生成するNumo::Randomを作った。アイディア自体はずっと前に考えていたが、腱鞘炎と忙しさで、形にできていなかった。

github.com

動機

Numo::NArray自体にも、一様乱数を生成する rand メソッドと、正規乱数を生成する rand_norm メソッドがある。さらに乱数のシードを指定する srand メソッドがある。シードを固定して、正規分布に従う乱数による大きさ100の配列を作るとすると次の様になる。

require 'numo/narray'

Numo::NArray.srand(42)
Numo::DFloat.new(100).rand_norm

とても簡単で便利だ。ただ、シードを固定する srand メソッドがグローバルに影響するので、あるメソッド内だけでシードを固定したいというのは難しい。

def hoge
  Numo::NArray.srand(42)
  Numo::DFloat.new(100).rand_norm
end

hoge

p Numo::DFloat.new(1).rand
# Numo::DFloat#shape=[1]
# [0.115514]

hoge

p Numo::DFloat.new(1).rand
# Numo::DFloat#shape=[1]
# [0.115514]

また、確率統計や機械学習のためには、一様分布や正規分布以外の確率分布も欲しくなるところ。

使い方

まず、Numo::Randomをインストールする。Numo::NArrayは依存Gemとして、一緒にインストールされる。

$ gem install numo-random

今回は、結果の可視化のために、Numo::Gnuplot(とgnuplot)をインストールする。

$ brew install gnuplot
$ gem install numo-gnuplot

Numo::Randomで、正規分布に従う乱数列の配列を作成する場合、次の様になる。

require 'numo/narray'
require 'numo/gnuplot'

require 'numo/random'

# 乱数生成器をシードを固定した作成する。
rng = Numo::Random::Generator.new(seed: 42)

# 平均0・分散1の正規分布に従う乱数による、大きさ (5000, 2) のNumo::DFloatを生成する。
x = rng.normal(shape: [5000, 2], loc: 0.0, scale: 1.0)

# 結果をプロットする。
Numo.gnuplot do
  set(terminal: 'png')
  set(output: 'normal2d.png')
  plot(x[true, 0], x[true, 1])
end

これを実行すると次の画像が得られる。中心が密で円形にサンプル点が広がっていて、正規分布だなという感じ。

Numo::Randomで作成した正規分布に従う乱数列

シードは乱数生成器のなかにあるので、グローバルな影響はない。

# シードを固定して乱数生成器を作成して、一様分布に従う乱数を得る。
rng_a = Numo::Random::Generator.new(seed: 8)
p rng_a.random
# 0.33322124834504996

# シードの異なる別の乱数生成器を生成しても、もとの乱数生成器が影響受けることはない。
rng_b = Numo::Random::Generator.new(seed: 42)
p rng_b.random
# 0.15802686859384152

rng_b = Numo::Random::Generator.new(seed: 42)
p rng_a.random
# 0.005744143787797302

一様分布や正規分布の他に、二項分布やガンマ分布に従う乱数生成メソッドを用意した。それはAPIドキュメントで。

yoshoku.github.io

実装

乱数生成アルゴリズムにはPermuted Congruential Generator (PCG) を利用した。高速で統計的性質が良いということで採用した。

www.pcg-random.org

PCGのC++実装が公式にあり、C++11から様々な確率分布に従う乱数生成ができるので、あとはこれらをnative extensionsでまとめ上げれば良いと考えた。

en.cppreference.com

Numo::NArrayな配列を生成した乱数で埋める共通ルーチンを書けば、あとは単純作業でメソッドを増やせる。少ないタイプ数で作れるので、腱鞘炎でも大丈夫だ。

おわりに

確率分布に従う乱数によるNumo::NArrayを生成するNumo::Randomを作った。今後のアップデートとしては、PCG以外の乱数生成アルゴリズムに対応することが考えられる。が、自分が欲しいと思っていたものはできたし、あとは手のダメージと相談といったところ。