Random numbers with vdRngGaussianMV.

[Fiori]

Hi!

I want to generate random numbers from a multivariate normal distribution with parameters mu = [3.0 5.0 2.0] and sigma = [ 16.0 8.0  4.0;  8.0 13.0 17.0; 4.0 17.0 62.0].

I wrote a mex-code using the the function vdRngGaussianMV and the results of the simulation are

>> mean(out)
ans =
    3.4558    2.9934    3.3013
>> cov(out)
ans =
  146.5081   -1.9068   -2.1787
   -1.9068  144.7461   -0.7771
   -2.1787   -0.7771  142.3955

Could you please tell me what I am doing wrong.

Thank you very much.

The code is the following:

#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)
#define BRNG    VSL_BRNG_MCG31 // VSL basic generator to be used

double normalMVN(int npars, int N, double *cov, double *mean, double *out);


/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    
    int npars,  N;
    double *cov, *mean, *out;
    
    /* make pointers to input data */
    npars = (int)mxGetScalar(prhs[0]);
    N = (int)mxGetScalar(prhs[1]);
    cov = mxGetPr(prhs[2]);
    mean = mxGetPr(prhs[3]);
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix( N, npars, mxREAL);
    out = mxGetPr(plhs[0]);
    
    
    /* call */
    normalMVN(npars, N, cov, mean, out);
    
    return 0;
    
}


double normalMVN(int npars, int N, double *cov, double *mean, double *out)
{
    
    /* This we need it for distributions */
    VSLSSTaskPtr task;
    VSLStreamStatePtr stream;
    
    char   uplo;
    int i, j, n, lda,  info;
    double *T, fiori[3];
    
    T = (double *)mxMalloc( npars*npars*sizeof( double* ) );
    
    uplo = 'L';
    n = npars;
    lda = npars;
    
    cblas_dcopy(npars*npars, cov, 1, T, 1);
    
    dpotrf( &uplo, &n, T, &lda, &info );
    if(info != 0){mexPrintf("c++ error: Cholesky failed\n\n");}
    
    vslNewStream( &stream, BRNG, SEED );
    
    vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2, stream, N, out, npars, VSL_MATRIX_STORAGE_FULL, mean, T );
  
    vslDeleteStream( &stream );
    
    /* Free memory */
    mxFree(T);
    
    return 0;
}

[Andrey_N_Intel]

Hi,

Per quick analysis of your code, it makes sense to apply a few corrections:

- replace uplo = 'L'; with uplo = 'U';

-remove '*' in sizeof( double* ) in T = (double *)mxMalloc( npars*npars*sizeof( double* ) );

For debugging purposes can you please set SEED  to a fixed number, provide value of N, and the data set you get from the generator for those parameters (I assume N is not huge, something like ~1000) ? It would help to compare the results on your and our side.

Also, which library version do you use?

Andrey

[Fiori]

Hello.

I corrected the code as you proposed.  I set the  SEED equal to 777,  npars =3, N=5000, sigma = [ 16.0 8.0  4.0;  8.0 13.0 17.0; 4.0 17.0 62.0];,  m=[3.0 5.0 2.0]; I have installed on my pc  "Intel Parallel Studio XE 2013 SP1, Update 2 for Windows".

I run the code with the command  out = fio(npars, N, sigma, m); and the results are

>> mean(out)
ans =
    3.3099    3.4070    3.3614
>> cov(out)
ans =
   32.5436   -0.5491   -0.9751
   -0.5491   29.8680   -0.6707
   -0.9751   -0.6707   32.0381
>>

Please see the attached txt for the simulated data.

The corrected code is:

#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 777
#define BRNG  VSL_BRNG_MCG31 // VSL basic generator to be used

double normalMVN(int npars, int N, double *cov, double *mean, double *out);


/* main fucntion */
void mexFunction(int nlhs,  mxArray *plhs[], int nrhs, const mxArray *prhs[])
{
    
    int npars,  N;
    double *cov, *mean, *out;
    
    /* make pointers to input data */
    npars = (int)mxGetScalar(prhs[0]);
    N = (int)mxGetScalar(prhs[1]);
    cov = mxGetPr(prhs[2]);
    mean = mxGetPr(prhs[3]);
    
    /* make pointers to output data */
    plhs[0] = mxCreateDoubleMatrix( N, npars, mxREAL);
    out = mxGetPr(plhs[0]);
    
    
    /* call */
    normalMVN(npars, N, cov, mean, out);
    
    return 0;
    
}


double normalMVN(int npars, int N, double *cov, double *mean, double *out)
{
    
    char   uplo;
    int i, j, n, lda,  info;
    double *T;
    
    
    /* This we need it for distributions */
    VSLSSTaskPtr task;
    VSLStreamStatePtr stream;
    vslNewStream( &stream, BRNG, SEED );
    
    
    /* Following variables are used in Cholesky factorization subroutine */
    T =  mxMalloc( npars*npars*sizeof( double ) );
    
    uplo = 'U';
    n = npars;
    lda = npars;
    
    cblas_dcopy(npars*npars, cov, 1, T, 1);
    
    /* MKL Choelsky factorization routine call */
    dpotrf( &uplo, &n, T, &lda, &info );
    if(info != 0){mexPrintf("c++ error: Cholesky failed\n\n");}
    
    
    /*****************************************************
     * Generating random numbers
     * from multivariate normal distribution
     *****************************************************/
    vdRngGaussianMV( VSL_METHOD_DGAUSSIANMV_BOXMULLER2, stream, N, out, npars, VSL_MATRIX_STORAGE_FULL, mean, T );
    
    
    vslDeleteStream( &stream );
    
    /* Free memory */
    mxFree(T);
    
    return 0;
}

[Andrey_N_Intel · Accepted solution]

Hi,

Thanks for the generated data. I reproduced covariance and mean estimates with your dataset. This indicates the bug in that version of Intel MKL. I attach C test case that was based on the your code, and the output produced by this test case. For the experiments I used Intel MKL from Intel Parallel Studio 2016 for Linux. Can you please download the latest version of the Studio and run the experiments with it?

Thanks, Andrey

[Fiori]

Hi!

I really appreciate your help. I installed on my pc   'Intel® Parallel Studio XE Cluster Edition , 2016 for Windows'  but  I don't know how to run my mex-code. With  "Intel Parallel Studio XE 2013 SP1" I used the command

'mex -v -g phi1c.c -I"C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl\include" -L"C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl\lib\intel64"  -lmkl_intel_lp64 -lmkl_core -lmkl_sequential -lmkl_cdft_core -lmkl_intel_thread'

Could you please tell me if you know  what to do?

Thank you very much.

[Ying_H_Intel]

Hi GF G., 

Are you asking the eqivalent folder in Intel® Parallel Studio XE Cluster Edition , 2016 for Windows? 

We changed the folder a little, for example, 

 C:\Program Files (x86)\Intel\Composer XE 2013 SP1\mkl

in 

C:\Program Files (x86)\Intel_sw_development_tools\parallel_studio_xe_2016.0.0.xx\compilers_and_libraries_2016\windows\mkl

You may replace them in your command 

Best Regards

Ying