2015-11-06 232 views
1

我想使用函数vdRngGaussian正常(平均值,sigma2)中绘制n个随机变量。 一种方式做,这是命令使用vdRngGaussian提出正态随机变量

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

取而代之的是我想用一个for循环。该MEX-代码我写的是

#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 

如果我把我的函数5次无一个for循环我得到这些结果

-0.2720, 2.1457, -1.2397, 0.7501, 0.1490 

你能帮我吗。 非常感谢。

回答

1

要创建与vslNewStream小溪种子(随机数发生器的初始状态)作为参数,而你就time function from the C time librarythe output基础此。问题是它返回。这个分辨率太粗糙了,当你在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过去的CPU时钟,但是当您关闭机器时,这会回到零。 或只是time in seconds as a double with MKL's dsecnd

+0

感谢您的帮助。我试图按照你的建议使用“high_resolution_clock”,但是当我运行代码时,我得到这个消息:“error:identifier”使用“未定义的名称空间标准;”。我试图添加库来修复它,但我没有运气。我读了函数“dscend”,并且它返回两次函数调用之间的时间(如果我理解正确的话)。你能否告诉我更多的细节如何使用它来定义一个新的种子?我很抱歉,也许这个问题很愚蠢,但我对C是新手。非常感谢。 –

+0

@ F.F。这意味着你不能编译为C++,或者你的编译器不支持C++ 11,或者你没有启用C++ 11支持。你能给编译器详细信息吗?您没有将问题标记为任何语言。对于'dsecnd',它通常用于计时,但它返回一个绝对时间,可能会被使用(乘以1e7左右,并舍入到要求的整数类型中)。 – chappjc

+0

你好@chappjc。我遇到了一些意想不到的问题,而且我以前也看不到yor的消息。我用扩展名.cpp保存了我的代码,它工作正常。非常感谢你。 –