采样变换 - rexp 与 rweibull

Question

我正在使用不同的采样函数，我想知道为什么这两种公式没有给出相同的结果

n=2
set.seed(1)
rweibull(n,shape = 1,scale = 1)
# [1] 1.3261078 0.9885284
set.seed(1)
rexp(n,rate = 1)
# [1] 0.7551818 1.1816428

当这等效时：

x <- c(0, rlnorm(50))
all.equal(dweibull(x, shape = 1), dexp(x))

是不是逆变换采样的问题？

如果是，为什么？

谢谢，

Answer 1

首先，

d*

表示分布函数，它是确定性的，这就是为什么你可以在下面的代码中获得相同的结果

x <- c(0, rlnorm(50))
all.equal(dweibull(x, shape = 1), dexp(x)

但是，

r*

给出给定分布的随机样本，这取决于随机性是如何触发的。

• rweibull

的源代码

#include "nmath.h"

double rweibull(double shape, double scale)
{
    if (!R_FINITE(shape) || !R_FINITE(scale) || shape <= 0. || scale <= 0.) {
    if(scale == 0.) return 0.;
    /* else */
    ML_ERR_return_NAN;
    }

    return scale * pow(-log(unif_rand()), 1.0 / shape);
}

• rexp

的源代码

#include "nmath.h"

double rexp(double scale)
{
    if (!R_FINITE(scale) || scale <= 0.0) {
    if(scale == 0.) return 0.;
    /* else */
    ML_ERR_return_NAN;
    }
    return scale * exp_rand(); // --> in ./sexp.c
}

其中

[exp_rand()][3]

用于生成指数随机变量，这与

pow(-log(unif_rand()), 1.0 / shape)

中所示的

rweibull.c

并不简单，而是更复杂，如下所示

#include "nmath.h"

double exp_rand(void)
{
    /* q[k-1] = sum(log(2)^k / k!)  k=1,..,n, */
    /* The highest n (here 16) is determined by q[n-1] = 1.0 */
    /* within standard precision */
    const static double q[] =
    {
    0.6931471805599453,
    0.9333736875190459,
    0.9888777961838675,
    0.9984959252914960,
    0.9998292811061389,
    0.9999833164100727,
    0.9999985691438767,
    0.9999998906925558,
    0.9999999924734159,
    0.9999999995283275,
    0.9999999999728814,
    0.9999999999985598,
    0.9999999999999289,
    0.9999999999999968,
    0.9999999999999999,
    1.0000000000000000
    };

    double a = 0.;
    double u = unif_rand();    /* precaution if u = 0 is ever returned */
    while(u <= 0. || u >= 1.) u = unif_rand();
    for (;;) {
    u += u;
    if (u > 1.)
        break;
    a += q[0];
    }
    u -= 1.;

    if (u <= q[0])
    return a + u;

    int i = 0;
    double ustar = unif_rand(), umin = ustar;
    do {
    ustar = unif_rand();
    if (umin > ustar)
        umin = ustar;
    i++;
    } while (u > q[i]);
    return a + umin * q[0];
}