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

29

如何在ruby中生成正态分布随机数的代码?

(注意:我回答了自己的问题,但我会等几天看看是否有更好的答案。)

编辑:

搜索此问题时,我查看了SO上所有与以下两个搜索结果相关的页面:

+"normal distribution" ruby

+gaussian +random ruby


你有查看相关问题吗?(请参见右侧面板) - Gumbo
是的,我检查过了,虽然有地方有这个算法,但没有人用Ruby编写过它。这是一个非常常见的任务,应该在标准库中。但如果没有的话,我认为可以在StackOverflow上找到复制粘贴的代码。 - Eponymous
最好提及你已经检查过哪些内容,这样想要回答的人就不会再次检查它们,除非他们认为你漏掉了某些东西。 - Andrew Grimm
4个回答

53

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 (CC0)

在法律允许的范围内,antonakos已经放弃了与RandomGaussian Ruby类相关的版权和邻近权利。此作品发表自:丹麦。


2
你能否给你的代码添加一个宽松的许可证(如BSD/CC-0或类似的)(因为它旨在供剪切粘贴重用)? - Eponymous
@antonakos 在我的OpenOffice Calc中,我有一个名为NORMDIST的函数,它有参数(Number, Mean, STDEV, C),我该如何使用您的代码来提供这4个参数?(您的代码接受'mean'、'stddev'和'rand',这与我在Excel/OpenOffice Calc中使用的不同) - medBouzid
我想知道为什么你要使用 Math.log(1 - rand.call) 而不是直接使用 Math.log(rand.call),因为 rand.call 返回的是 0 到 1 之间的数字。我已经检查了 Python 和 Boost 代码,它们也是这样做的...但是为什么呢? - kaikuchn
Ruby中最小的非零浮点数是0.0.next_float => 5.0e-324,小于1.0的最大浮点数是1.0.prev_float => 0.9999999999999999。我不知道Kernel.rand是否会返回5.0e-324,但如果它确实这样做,那么减法将截断随机变量的尾部。 - kaikuchn

21
原问题要求提供代码,但作者的后续评论暗示其对使用现有库也很感兴趣。我也对此感兴趣,并进行了搜索,找到了以下两个Ruby宝石: gsl -“GNU科学图书馆的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

12

+1 持 @antonakos 的观点。以下是我一直在使用的 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 Algorithm:)。


12

另一个选项是使用distribution gem,这个gem是由SciRuby fellows之一编写的。

我认为它使用起来更简单。

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

看这段代码,我完全不知道预期标准差应该是多少。 - Paul Brannan
2
请访问此处:http://rubydoc.info/gems/distribution/0.7.0/Distribution/Normal/Ruby_您会发现rng(1)指定了平均值,并且您可以通过传递附加参数Distribution::Normal.rng(mean, standard_deviation)来指定所需的标准差,然后调用它以提供该分布中的随机值。 - Ryanmt

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接