[转载]从多维正态分布中随机抽取样本

在PCL的VoxelGridCovariance类的getDisplayCloud方法中采用了Cholesky分解采样的方法。

template<typename PointT> void
pcl::VoxelGridCovariance<PointT>::getDisplayCloud (pcl::PointCloud<PointXYZ>& cell_cloud)
{
  cell_cloud.clear ();

  int pnt_per_cell = 1000;
  boost::mt19937 rng;
  boost::normal_distribution<> nd (0.0, leaf_size_.head (3).norm ());
  boost::variate_generator<boost::mt19937&, boost::normal_distribution<> > var_nor (rng, nd);

  Eigen::LLT<Eigen::Matrix3d> llt_of_cov;
  Eigen::Matrix3d cholesky_decomp;
  Eigen::Vector3d cell_mean;
  Eigen::Vector3d rand_point;
  Eigen::Vector3d dist_point;

  // Generate points for each occupied voxel with sufficient points.
  for (typename std::map<size_t, Leaf>::iterator it = leaves_.begin (); it != leaves_.end (); ++it)
  {
    Leaf& leaf = it->second;

    if (leaf.nr_points >= min_points_per_voxel_)
    {
      cell_mean = leaf.mean_;
      llt_of_cov.compute (leaf.cov_);
      cholesky_decomp = llt_of_cov.matrixL ();

      // Random points generated by sampling the normal distribution given by voxel mean and covariance matrix
      for (int i = 0; i < pnt_per_cell; i++)
      {
        rand_point = Eigen::Vector3d (var_nor (), var_nor (), var_nor ());
        dist_point = cell_mean + cholesky_decomp * rand_point;
        cell_cloud.push_back (PointXYZ (static_cast<float> (dist_point (0)), static_cast<float> (dist_point (1)), static_cast<float> (dist_point (2))));
      }
    }
  }
}

---

原文链接： http://blog.sina.com.cn/s/blog_955cedd8010130m8.html

 

R = mvnrnd(MU,SIGMA)——从均值为MU，协方差为SIGMA的正态分布中抽取n*d的矩阵R（n代表抽取的个数，d代表分布的维数）。

MU为n*d的矩阵，R中的每一行为以MU中对应的行为均值的正态分布中抽取的一个样本。

SIGMA为d*d的对称半正定矩阵，或者为d*d*n的array。若SIGMA为array，R中的每一行对应的分布的协方差矩阵为该array对应的一个page。也就是说：R(i,:)由MU(i,:)和SIGMA(:,:,i)产生。
如果协方差矩阵为对角阵，sigma也可用1*d向量或1*d*n的array表示，如果MU是一个1*d的向量，则SIGMA中的n个协方差矩阵共用这个MU。R的行数n由MU的行数n或者SIGMA的page数n决定。

r = mvnrnd(MU,SIGMA,cases)——从均值为MU(1*d)，协方差矩阵为SIGMA(d*d)的正态分布中随机抽取cases个样本，返回cases*d的矩阵r。

 

不使用现成的函数，可以通过一个线性变换来实现：

我们知道，matlab产生的n维正态样本中的每个分量都是相互独立的，或者说，它的协方差矩阵是一个数量矩阵mI，如：X = randn(10000,4);产生10000个4维分布的正态分布样本，协方差矩阵就是单位矩阵I。

定理 n维随机变量X服从正态分布N(u,B)，若m维随机变量Y是X的线性变换，即Y=XC，其中C是n×m阶矩阵，则Y服从m维正态分布N(uC,C'BC)。

根据这条定理，我们可以通过一个线性变换C把协方差矩阵为I的n维正态样本变为协方差矩阵为C'C的n维正态样本。如果我们要产生协方差矩阵为R的n维正态样本，由于R为对称正定矩阵，所以有Cholesky分解: R=C'C

附：matlab程序

function Y = multivrandn(u,R,M)
% this function draws M samples from N(u,R)
% where u is the mean vector(row) and R is the covariance matrix which must be positive definite

n = length(u);              % get the dimension
C = chol(R);                % perform cholesky decomp R = C'C
X = randn(M,n);             % draw M samples from N(0,I)

Y = X*C + ones(M,1)*u;