抽样变换REXP与RWEIBUR

发布于01月31日

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

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

是逆变换抽样的问题吗？

若有，原因为何？

谢谢,

推荐答案

首先，d*表示密度函数，它是确定性的，这就是为什么您可以在以下代码中获得相同的结果

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

但是，r*会提供给定分布的随机样本，这取决于随机性是如何触发的.

百人组

#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);
}

百人组

#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()是用来产生指数随机变量的，它不是简单地与rweibull.c中的pow(-log(unif_rand()), 1.0 / shape)相同，而是更复杂，如下所示

#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];
}