将均匀分布转换为正态分布

130

如何将均匀分布(大多数随机数生成器产生的分布,例如在0.0和1.0之间)转换为正态分布?如果我想要选择自己的均值和标准差怎么办?


3
你有语言规范吗,还是这只是一个通用的算法问题? - Bill the Lizard
3
一般算法问题。我不在乎使用哪种语言,但我希望答案不依赖于特定语言提供的特定功能。 - Terhorst
16个回答

58

有很多方法:

  • 不要使用Box Muller。尤其是当你需要产生许多高斯数时。Box Muller生成的结果被夹在-6到6之间(假设采用双精度浮点数,使用单精度浮点数效果更差)。而且它比其他可用的方法要低效得多。
  • Ziggurat可以使用,但需要进行表查找(由于缓存大小问题,需要进行一些特定于平台的调整)。
  • Ratio-of-uniforms是我最喜欢的方法,只需要进行一些加法/乘法运算,并且1/50的时间内进行对数运算(例如,请参考此处)。
  • 反转CDF确实是有效的(为什么被忽视了呢?),如果你搜索谷歌,可以找到快速实现。对于准随机数来说,这是必须的。

3
你确定要使用[-6,6]夹紧吗?如果属实,这是一个非常重要的点(值得在维基百科页面上注明)。 - redcalx
9
@locster:这种不良特性也存在于逆CDF方法中。请参阅http://www.cimat.mx/~src/prope08/randomgauss.pdf。可以通过使用具有接近零的浮点数的非零概率的均匀RNG来缓解此问题。大多数RNG不会这样做,因为它们生成一个(通常为64位)整数,然后将其映射到[0,1]。这使得这些方法不适用于对高斯变量的尾部进行采样(例如在计算金融中定价低/高行权期权)。 - Alexandre C.
1
考虑 x1 = 1 / (2^32-1)(使用32位整数生成的大于零的最小随机数),x2 = 0。然后,random = sqrt(-2 ln(1/(2^32-1))) * cos(2 * pi * 0) = 6.6604。使用64位随机数生成器,这个值约为9.419。 - the swine
11
@AlexandreC。为了明确两点,使用64位数字,尾巴会延伸到8.57或9.41(将其转换为[0,1)后取对数的较低值)。即使夹在[-6, 6]范围内,超出此范围的机会约为1.98e-9,对于大多数人来说已经足够好了,即使是在科学领域。对于8.57和9.41这些数字,这个概率变成了1.04e-17和4.97e-21。这些数字非常小,以至于在上述限制条件下,使用Box-Muller采样和真正的高斯采样之间的差别几乎纯粹是学术性的。如果需要更好的效果,只需将它们累加四次并除以2即可。 - CrazyCasta
9
我认为建议不使用Box Muller变换会误导大多数用户。了解其限制是很好的,但正如CrazyCasta所指出的,对于大多数不严重依赖于异常值的应用程序,你可能不需要担心这个问题。例如,如果你曾经依赖于从numpy中的正态分布进行抽样,那么你就依赖于Box Muller变换(极坐标形式)。https://github.com/numpy/numpy/blob/c08d2647240555e730da7580374a61d8547a932e/numpy/random/mtrand/randomkit.c#L619。 - Andreas Grivas
显示剩余7条评论

51

7
这些方法都适用于线性同余生成器的常见警告,因此请使用良好的底层生成器。干杯。 - dmckee --- ex-moderator kitten
3
像Mersenne Twister这样的算法,还是您有其他建议? - Gregg Lind

34

将任何函数的分布转换为另一个函数的分布,需要使用您想要的函数的反函数。

换句话说,如果您想要特定的概率函数 p(x),则通过对其积分-> d(x) = integral(p(x))并使用其反函数:Inv(d(x))来获得分布。现在使用随机概率函数(具有均匀分布),并将结果值通过函数 Inv(d(x)) 转换。您应该得到根据您选择的函数分布投掷的随机值。

这是一种通用的数学方法 - 使用它,您现在可以选择您拥有的任何概率或分布函数,只要其具有反函数或良好的反函数逼近。

希望这有所帮助,并感谢关于使用分布而不是概率本身的小提醒。


4
这是一种被忽视的生成高斯变量的方法,效果非常好。在这种情况下,可以使用牛顿法有效地计算反向累积分布函数(导数为e^{-t^2}),并且很容易得到一个有理数近似值作为初始估计,因此您需要进行3-4次erf和exp的评估。如果您使用准随机数,则必须使用正好一个均匀数来获取高斯数,这是必需的。 - Alexandre C.
10
请注意,你需要反转累积分布函数而不是概率密度函数。Alexandre已经暗示了这一点,但我认为更明确地提到这一点可能有助于理解——因为答案似乎表明了PDF。 - ltjax
如果你准备相对于平均值随机选择一个方向,那么你可以使用PDF吗?我理解得对吗? - Mark McKenna
2
这被称为反变换抽样 - dashesy
2
这里是SE上一个相关问题,有更通用的答案和很好的解释。 - dashesy

23

这是使用 Box-Muller 转换的极坐标形式的 JavaScript 实现。

/*
 * Returns member of set with a given mean and standard deviation
 * mean: mean
 * standard deviation: std_dev 
 */
function createMemberInNormalDistribution(mean,std_dev){
    return mean + (gaussRandom()*std_dev);
}

/*
 * Returns random number in normal distribution centering on 0.
 * ~95% of numbers returned should fall between -2 and 2
 * ie within two standard deviations
 */
function gaussRandom() {
    var u = 2*Math.random()-1;
    var v = 2*Math.random()-1;
    var r = u*u + v*v;
    /*if outside interval [0,1] start over*/
    if(r == 0 || r >= 1) return gaussRandom();

    var c = Math.sqrt(-2*Math.log(r)/r);
    return u*c;
    
    /* todo: optimize this algorithm by caching (v*c) 
     * and returning next time gaussRandom() is called.
     * left out for simplicity */
}

你的r变量的区间实际上是(0,1],而不是[0,1],因为当r等于零时,你会重新开始。 - jpaugh

6

其中R1,R2是随机均匀数:

正态分布,标准差为1:

sqrt(-2*log(R1))*cos(2*pi*R2)

这就是精确的方法...不需要进行那些缓慢的循环!

参考资料:dspguide.com/ch2/6.htm


在有人纠正之前...这是我想出来的近似值:(1.5 - (R1+R2+R3)) * 1.88。我也喜欢它。 - Erik Aronesty
1
谢谢,我也在这里找到了这个方程式 http://www.dspguide.com/ch2/6.htm - iperov
2
这是其他答案中提到的Box-Muller变换,具有相同的限制条件。请注意,您可以通过计算正弦值来从R1和R2获得第二个正态随机偏差。有关Box-Muller的更多详细信息,请参见此处:https://www.baeldung.com/cs/uniform-to-normal-distribution。 - ELNJ

6

1
这不会产生特别接近的正态分布(“尾部”或端点不会接近真实的正态分布)。Box-Muller方法更好,就像其他人建议的那样。 - Peter K.
1
Box Muller 的尾部也是错误的(它以双精度返回-6到6之间的数字)。 - Alexandre C.
1
n=12(在0到1的范围内求和12个随机数,并减去6)的结果为stddev=1和mean=0。然后可以使用它来生成任何正态分布。只需将结果乘以所需的stddev并加上mean即可。 - JerryM

4

八年后我仍能为Java添加一些内容,特别是关于Random.nextGaussian()方法,该方法可为您生成均值为0.0、标准差为1.0的高斯分布。

只需进行简单的加法和乘法运算即可将其平均值和标准差更改为您所需的数值。


4
我会使用Box-Muller方法。需要注意以下两点:
  1. 每次迭代会产生两个值
    通常,您将一个值存储在缓存中,并返回另一个值。在下一次调用时,您返回缓存的值。
  2. Box-Muller方法会生成Z-score值
    你需要通过标准差来缩放Z-score,并加上均值,才能得到正态分布中的完整值。

你如何扩展Z分数? - Terhorst
3
缩放 = 平均值 + 标准差 * z分数 // 得到正态分布(mean, stdDev^2) - yoyoyoyosef

1
这是一个使用Box-Muller变换的极坐标形式的Matlab实现:
函数randn_box_muller.m:
function [values] = randn_box_muller(n, mean, std_dev)
    if nargin == 1
       mean = 0;
       std_dev = 1;
    end

    r = gaussRandomN(n);
    values = r.*std_dev - mean;
end

function [values] = gaussRandomN(n)
    [u, v, r] = gaussRandomNValid(n);

    c = sqrt(-2*log(r)./r);
    values = u.*c;
end

function [u, v, r] = gaussRandomNValid(n)
    r = zeros(n, 1);
    u = zeros(n, 1);
    v = zeros(n, 1);

    filter = r==0 | r>=1;

    % if outside interval [0,1] start over
    while n ~= 0
        u(filter) = 2*rand(n, 1)-1;
        v(filter) = 2*rand(n, 1)-1;
        r(filter) = u(filter).*u(filter) + v(filter).*v(filter);

        filter = r==0 | r>=1;
        n = size(r(filter),1);
    end
end

调用histfit(randn_box_muller(10000000),100);的结果如下图所示: Box-Muller Matlab Histfit

很明显,与Matlab内置函数randn相比,这种方法效率非常低。


1

标准的Python库模块random拥有您所需的功能:

normalvariate(mu, sigma)
正态分布。mu是平均值,sigma是标准差。

关于算法本身,请查看Python库中random.py中的函数。

手册条目在此处


2
不幸的是,Python的库使用了Kinderman,A.J.和Monahan,J.F.的算法,“Computer generation of random variables using the ratio of uniform deviates”,ACM Trans Math Software, 3, (1977), pp257-260。该算法使用两个均匀分布的随机变量来生成正态值,而不是一个,因此不太明显如何将其用作OP所需的映射。 - Ian

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