关于算法:将均匀分布转换为正态分布

Converting a Uniform Distribution to a Normal Distribution

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


Ziggurat算法对此非常有效,尽管Box-Muller变换更容易从头开始实现(而不是疯狂慢)。


有很多方法:

  • 请勿使用Box Muller。特别是当您绘制许多高斯数时。 Box Muller得出的结果被限制在-6和6之间(假设双精度。浮点会使情况更糟)。而且它确实比其他可用方法效率低。
  • Ziggurat很好,但需要进行表查找(由于缓存大小问题,需要进行一些特定于平台的调整)
  • 均匀率是我的最爱,只有几次加法/乘法和对数的1/50(例如,看那里)。
  • 反转CDF是高效的(并且被忽略了,为什么?),如果您搜索google,就可以快速实现它。准随机数是必需的。


将任何函数的分布更改为另一个函数都涉及使用所需函数的逆函数。

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

这是通用的数学方法-通过使用它,您现在可以选择具有逆或良好逆近似的任何概率或分布函数。

希望这对您有所帮助,并感谢您对使用分布的简短评论,而不是概率本身。


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

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
/*
 * 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 */
}

使用中央极限定理Wikipedia条目mathworld条目可发挥您的优势。

生成n个均匀分布的数字,将它们相加,减去n * 0.5,您将得到近似正态分布的输出,均值等于0,方差等于(1/12) * (1/sqrt(N))(关于最后一个均匀分布,请参见Wikipedia)

n = 10可以使您快得快一半。如果您想要的东西超过一半,请寻求轮胎解决方案(正态分布在Wikipedia条目中已指出)


我会使用Box-Muller。关于以下两点:

  • 您每次迭代都会得到两个值
    通常,您缓存一个值,然后返回另一个值。在下一次调用样本时,您返回缓存的值。
  • Box-Muller给出Z得分
    然后,您必须按标准偏差缩放Z分数,并添加均值以获得正态分布中的完整值。

  • 其中R1,R2是随机统一数字:

    正态分布,SD为1:sqrt(-2 * log(R1))* cos(2 * pi * R2)

    这是正确的……不需要执行所有这些慢循环!


    八年后我可以添加一些东西似乎令人难以置信,但是对于Java,我想向读者介绍Random.nextGaussian()方法,该方法为您生成均值0.0和标准差1.0的高斯分布。

    简单的加法和/或乘法将改变均值和标准差,以满足您的需求。


    标准的Python库模块random具有您想要的:

    normalvariate(mu, sigma)
    Normal distribution. mu is the mean, and sigma is the standard deviation.

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

    手动输入在这里


    问:如何将均匀分布(如大多数随机数生成器产生的,例如0.0到1.0之间)转换为正态分布?

  • 对于软件实现,我知道几个随机生成器名称,这些名称在[0,1]中为您提供了伪统一的随机序列(Mersenne Twister,线性一致生成器)。我们称它为U(x)

  • 存在于数学领域中称为概率论。
    第一件事:如果您想对r.v建模积分分布为F的情况下,您可以尝试仅求F ^ -1(U(x))的值。在理论上证明了这样的r.v.将具有积分分布F。

  • 当F ^ -1可以无问题地解析得出时,第2步可适用于生成r.v.?F,而无需使用任何计数方法。 (例如exp.distribution)

  • 为了建模正态分布,您可以计算y1 * cos(y2),其中y1?[2pi]是均匀的。 y2是相对分布。

  • 问:如果我要选择平均值和标准差怎么办?

    您可以计算sigma * N(0,1)+ m。

    可以证明,这种移位和缩放导致N(m,sigma)


    我想您应该在EXCEL:=norminv(rand();0;1)中尝试此操作。这将产生应该以零均值和单位方差正态分布的随机数。可以提供任何值" 0",以便数字具有期望的均值,并且通过更改" 1",您将获得与输入平方相等的方差。

    例如:=norminv(rand();50;3)将产生为MEAN = 50 VARIANCE = 9的正态分布数字。


    使用实现的函数rnorm()也比创建正态分布的随机数生成器要快,因为它比写随机数生成器要快。参见以下代码作为证明

    1
    2
    3
    4
    5
    n <- length(z)
    t0 <- Sys.time()
    z <- rnorm(n)
    t1 <- Sys.time()
    t1-t0

    这是使用Box-Muller变换的极坐标形式的Matlab实现:

    功能randn_box_muller.m

    1
    2
    3
    4
    5
    6
    7
    8
    9
    10
    11
    12
    13
    14
    15
    16
    17
    18
    19
    20
    21
    22
    23
    24
    25
    26
    27
    28
    29
    30
    31
    32
    33
    34
    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
    2
    3
    4
    5
    6
    7
    8
    set.seed(123)
    n <- 1000
    u <- runif(n) #creates U
    x <- -log(u)
    y <- runif(n, max=u*sqrt((2*exp(1))/pi)) #create Y
    z <- ifelse (y < dnorm(x)/2, -x, NA)
    z <- ifelse ((y > dnorm(x)/2) & (y < dnorm(x)), x, z)
    z <- z[!is.na(z)]

    1
    2
    3
    4
    5
    6
    function distRandom(){
      do{
        x=random(DISTRIBUTION_DOMAIN);
      }while(random(DISTRIBUTION_RANGE)>=distributionFunction(x));
      return x;
    }

    近似:

    1
    2
    3
    function rnd_snd() {
        return (Math.random()*2-1)+(Math.random()*2-1)+(Math.random()*2-1);
    }

    参见http://www.protonfish.com/random.shtml


    推荐阅读