Crystalで正規分布に基づく乱数を生成する

はじめに

Crystal言語で正規乱数が出したくなった。しかし、Crystalは数値計算を目的としたプログラミング言語ではないので標準ライブラリにはなかった。

Box-Muller法

ボックス=ミュラー法というものを使えば、一様分布に従う確率変数から、標準正規分布に従う確率変数を生成できるそうだ。

実装

Seedを設定して毎回同じ乱数が出るようにしたかった。

• 標準ライブラリのRandomをオープンクラスで拡張する
• 標準ライブラリのRandomを継承してRandom2を作る
• 新しくRandクラスを作り、知らないメソッドが呼ばれたらRandomに委譲する

の3つのパターンが考えられて、一番最後のものを採用することにした。

class Rand
  def initialize
    r = Random.new
    initialize(random: r)
  end

  def initialize(seed : UInt64)
    r = Random.new(seed)
    initialize(random: r)
  end

  def initialize(random : Random)
    @random = random
    @iset = 0.to_i32
    @gset = 0.to_f64
  end

  forward_missing_to @random

  # Generates a random number from a standard normal distribution.
  def randn : Float64
    if @iset == 0
      rsq, v1, v2, fac = 0.0, 0.0, 0.0, 0.0
      loop do
        v1 = 2.0 * @random.next_float - 1.0
        v2 = 2.0 * @random.next_float - 1.0
        rsq = v1 * v1 + v2 * v2
        break if rsq < 1.0 && rsq != 0.0
      end

      fac = Math.sqrt(-2.0 * Math.log(rsq) / rsq)
      @gset = v1 * fac
      @iset = 1
      return v2 * fac
    else
      @iset = 0
      return @gset
    end
  end

  # Generates a random number from a normal distribution with the specified mean and standard deviation.
  def randn(mean, std_dev) : Float64
    mean = mean.to_f
    std_dev = std_dev.to_f
    raise "std_dev must be positive" if std_dev.negative?
    mean + std_dev * randn()
  end
end

Cのコードなどを見ながら見よう見まねで書いている。何か間違っている部分を発見したらコメント欄で指摘してほしい。

ライブラリにした

このコードは再利用する機会がありそうなのでShardにした。Crystalは発展途上の言語なので、自分で書いたコードを再利用しやすいように積極的にライブラリとして公開していきたい。

この記事は以上です。