可以通过GNU Scientific Library (GSL)实现。函数为gsl_ran_multivariate_gaussian ()
,具体细节可以参考官方文档
举个例子
$$ \mu = [1, 2, 3]\qquad \Sigma = diag\{1, 2, 3\} $$
可以通过下列C++程序实现
#include<iostream>
#include<gsl/gsl_vector.h>
#include<gsl/gsl_matrix.h>
#include<gsl/gsl_linalg.h>
#include<gsl/gsl_rng.h>
#include<gsl/gsl_randist.h>
#include<gsl/gsl_cdf.h>
using namespace std;
int main()
{
// setup mean and variance-covariance matrix
double mu[] = {1, 2, 3};
double omega[] = {1, 0, 0, 0, 2, 0, 0, 0, 3};
gsl_matrix_view v_omega = gsl_matrix_view_array(omega, 3, 3);
gsl_vector_view v_mu = gsl_vector_view_array(mu, 3);
// setup random generator
const gsl_rng_type *T;
gsl_rng *r;
gsl_rng_env_setup();
T = gsl_rng_default;
r = gsl_rng_alloc(T);
// cholesky decomp
gsl_linalg_cholesky_decomp(&v_omega.matrix);
// mvrnorm
gsl_vector *res = gsl_vector_calloc(3);
gsl_ran_multivariate_gaussian(r, &v_mu.vector, &v_omega.matrix, res);
for (int i = 0; i < 3; i++)
cout << gsl_vector_get(res, i) << " ";
cout << endl;
gsl_rng_free(r);
return 0;
}
结果为