box_muller.rb

require "gnuplot"

def draw_chart(x, y)
  Gnuplot.open do |gp|
    Gnuplot::Plot.new(gp) do |plot|
      y.each do |name, value|
        if x.size == value[0].size
          plot.data << Gnuplot::DataSet.new([x, value[0]]) do |ds|
            ds.with = value[1]
            ds.title = name
          end
        end
      end
    end
  end
end

def make_hist(x_array)
  bin = 1000

  min = x_array.min
  max = x_array.max
  len = max - min
  width = len / bin

  hist = Array.new(bin, 0)
  x_array.each_with_index do |x, i|
    # 最小値を0にしてからビンの位置を計算
    pos = ((x + min.abs) / width).to_i
    pos = bin - 1 if pos == bin
    hist[pos] += 1
  end

  sum = hist.inject(0.0){|s, h| s += h * width}
  hist.map! do |h|
    h /= sum
  end

  # 最小値を戻す
  chart_x = []
  bin.times do |i|
    chart_x << i * width - min.abs
  end

  # 正規分布
  normal_y = []
  chart_x.each do |x|
    normal_y << 1 / Math.sqrt(2 * Math::PI) * Math.exp(-x ** 2 / 2)
  end

  draw_chart(chart_x, {"Normal Distribution" => [normal_y, "lines"],
                       "Box Muller" => [hist, "points"]})
end

gaussian_variables = []
1E7.to_i.times do |i|
  gaussian_variables << Math.sqrt(-2 * Math.log(rand)) * Math.cos(2 * Math::PI * rand)
end
make_hist gaussian_variables