用于在Ruby中生成高斯(正态分布)随机数的代码

Question

什么是在ruby中生成正态分布随机数的代码？

(注意:我回答了我自己的问题,但我会等几天才接受是否有人有更好的答案.)

编辑:

搜索到这一点,我查看了两个搜索产生的SO上的所有页面:

+"正态分布"红宝石

和

+高斯+随机红宝石

Answer 1

Python的random.gauss()和Boost的normal_distribution都使用Box-Muller变换,因此对Ruby也应该足够好.

def gaussian(mean, stddev, rand)
  theta = 2 * Math::PI * rand.call
  rho = Math.sqrt(-2 * Math.log(1 - rand.call))
  scale = stddev * rho
  x = mean + scale * Math.cos(theta)
  y = mean + scale * Math.sin(theta)
  return x, y
end

该方法可以包含在一个逐个返回样本的类中.

class RandomGaussian
  def initialize(mean, stddev, rand_helper = lambda { Kernel.rand })
    @rand_helper = rand_helper
    @mean = mean
    @stddev = stddev
    @valid = false
    @next = 0
  end

  def rand
    if @valid then
      @valid = false
      return @next
    else
      @valid = true
      x, y = self.class.gaussian(@mean, @stddev, @rand_helper)
      @next = y
      return x
    end
  end

  private
  def self.gaussian(mean, stddev, rand)
    theta = 2 * Math::PI * rand.call
    rho = Math.sqrt(-2 * Math.log(1 - rand.call))
    scale = stddev * rho
    x = mean + scale * Math.cos(theta)
    y = mean + scale * Math.sin(theta)
    return x, y
  end
end

(CC0)

在法律允许的范围内,antonakos已放弃RandomGaussianRuby类的所有版权及相关或邻接权利.这项工作发表于:丹麦.

许可声明并不意味着我关心这段代码.相反,我不使用代码,我没有测试它,我不用Ruby编程.

Answer 2

原始问题要求代码,但作者的后续评论意味着对使用现有库感兴趣.我对此感兴趣,我的搜索出现了这两个红宝石宝石:

gsl - "GNU Scientific Library的Ruby接口"(要求您安装GSL).具有均值= 0和给定标准偏差的正态分布随机数的调用序列是

 rng = GSL::Rng.alloc
 rng.gaussian(sd)      # a single random sample
 rng.gaussian(sd, 100) # 100 random samples

rubystats - "来自PHPMath的统计库的一个端口"(纯ruby).具有给定均值和标准差的正态分布随机数的调用序列是

 gen = Rubystats::NormalDistribution.new(mean, sd)
 gen.rng               # a single random sample
 gen.rng(100)          # 100 random samples

Answer 3

@ antonakos答案+1.这是我一直在使用的Box-Muller的实现; 它本质上是相同但稍微紧凑的代码:

class RandomGaussian
  def initialize(mean = 0.0, sd = 1.0, rng = lambda { Kernel.rand })
    @mean, @sd, @rng = mean, sd, rng
    @compute_next_pair = false
  end

  def rand
    if (@compute_next_pair = !@compute_next_pair)
      # Compute a pair of random values with normal distribution.
      # See http://en.wikipedia.org/wiki/Box-Muller_transform
      theta = 2 * Math::PI * @rng.call
      scale = @sd * Math.sqrt(-2 * Math.log(1 - @rng.call))
      @g1 = @mean + scale * Math.sin(theta)
      @g0 = @mean + scale * Math.cos(theta)
    else
      @g1
    end
  end
end

当然,如果你真的关心速度,你应该实现Ziggurat算法 :).

Answer 4

另一个选项,这个使用分发宝石,由一个SciRuby研究员写的.

我认为它使用起来有点简单.

require 'distribution'
normal = Distribution::Normal.rng(1)
norm_distribution = 1_000.times.map {normal.call}