使用 vdRngGaussian 提出正态随机变量

Question

我想使用函数 vdRngGaussian 从 Normal(mean, sigma2) 中绘制 n 个随机变量。一种方法是使用命令

vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, stream, n, x, mean, sqrt(sigma2) )

我想使用 for 循环来代替这个。我写的 mex-code 是

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include <mkl.h>
#include "mkl_vml.h"
#include "mex.h"
#include "matrix.h"
#include "mkl_vsl.h"
#include <time.h>
#define SEED  time(NULL)

double normal(double mean, double sigma2, VSLStreamStatePtr stream);

/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    double *x, mean, sigma2;

    VSLStreamStatePtr stream;
    vslNewStream( &stream, VSL_BRNG_MT19937, SEED );

    /* make pointers to input data */
    mean = (double)mxGetScalar(prhs[0]);
    sigma2= (double)mxGetScalar(prhs[1]);

    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix(1, 1, mxREAL);
    x = mxGetPr(plhs[0]);


    x[0] = normal( mean, sigma2, stream);


    /* Deleting the stream */
    vslDeleteStream( &stream );

    return;
}

double normal(double mean, double sigma2, VSLStreamStatePtr stream)
{
    double x[1];

    vdRngGaussian(VSL_RNG_METHOD_GAUSSIAN_ICDF, stream, 1, x, mean, sqrt(sigma2) );

    return(x[0]);
}

当我运行代码用命令

for i=1:5
out(i) = normalc(0.0, 1.0);
end

我得到以下结果：

-1.1739   -1.1739   -1.1739   -1.1739   -1.1739

如果我在没有 for 循环的情况下调用我的函数 5 次，我会得到这些结果

-0.2720, 2.1457, -1.2397, 0.7501, 0.1490

你能帮帮我吗？非常感谢。

Answer 1

您正在使用 vslNewStream takes a seed (the initial condition of the random number generator) as an argument, and you are basing this on the output of the time function from the C time library 创建的流。问题是 returns 秒。这个分辨率太粗糙了，当你在 for 循环中调用它时，你最终会得到相同的种子。（假设 for 循环执行得很快，而且你很幸运，开始时间。）

也许使用不同的种子，例如 high_resolution_clock from the C++11 <chrono> header。您可能仍会得到纪元时间，但单位不同：

using namespace std::chrono;
system_clock::time_point now = high_resolution_clock::now();
system_clock::duration epoch = now.time_since_epoch();

long long ns = duration_cast<nanoseconds>(epoch).count();
long long mic = duration_cast<microseconds>(epoch).count();
long long ms = duration_cast<milliseconds>(epoch).count();

Demo (cpp.sh).

或者，您可以使用来自 MKL's mkl_get_cpu_clocks function, but this goes back to zero when you shutdown your machine. Or just the time in seconds as a double with MKL's dsecnd 的已用 CPU 时钟。

使用 vdRngGaussian 提出正态随机变量

Propose normal random variables using vdRngGaussian

mex

intel-mkl