产生满足二维高斯（正态）分布的随机数：

mu=[0,2];%数学期望
sigma=[1 0;0,4];%协方差矩阵
r=mvnrnd(mu,sigma,50)%生成50个样本

help mvnrnd

MVNRND Random vectors from the multivariate normal distribution.
    R = MVNRND(MU,SIGMA) returns an N-by-D matrix R of random vectors
    chosen from the multivariate normal distribution with mean vector MU,
    and covariance matrix SIGMA.  MU is an N-by-D matrix, and MVNRND
    generates each row of R using the corresponding row of MU.  SIGMA is a
    D-by-D symmetric positive semi-definite matrix, or a D-by-D-by-N array.
    If SIGMA is an array, MVNRND generates each row of R using the
    corresponding page of SIGMA, i.e., MVNRND computes R(I,:) using MU(I,:)
    and SIGMA(:,:,I).  If the covariance matrix is diagonal, containing
    variances along the diagonal and zero covariances off the diagonal,
    SIGMA may also be specified as a 1-by-D matrix or a 1-by-D-by-N array,
    containing just the diagonal. If MU is a 1-by-D vector, MVNRND
    replicates it to match the trailing dimension of SIGMA.
 
    R = MVNRND(MU,SIGMA,N) returns a N-by-D matrix R of random vectors
    chosen from the multivariate normal distribution with 1-by-D mean
    vector MU, and D-by-D covariance matrix SIGMA.
 
    Example:
 
       mu = [1 -1]; Sigma = [.9 .4; .4 .3];
       r = mvnrnd(mu, Sigma, 500);
       plot(r(:,1),r(:,2),'.');

生成二维高斯（正态）分布图形——meshgrid

%two-dimensional Normal Distribution
%
% (C)2008 TangSheng
%
function NormDis(u1,u2,sig1,sig2,rho)

�fault
if nargin<5, rho = 0; end                   %相关系数

if nargin<4, sig2 = 2;end                   %正态分布2方差
if nargin<3, sig1 = 1;end                   %正态分布1方差

if nargin<2, u2 = 2;end                     %正态分布2均值
if nargin<1, u1 = 6;end                     %正态分布1均值
 
%global ava;
ava = [u1,u2];                                %高斯分布均值向量
 
cov_xy = rho*sig1*sig2;                       %协方差
 
%global sigma;
sigma = [sig1 cov_xy ;cov_xy sig2 ];          %协方差矩阵
 
%------数据显示网格范围------------%
 
scop1_l = u1-sqrt(sig1)-2;
scop1_r = u1+sqrt(sig1)+2;
scop2_l = u2-sqrt(sig2)-2;
scop2_r = u2+sqrt(sig2)+2;
[X,Y] = meshgrid(scop1_l:0.2:scop1_r,scop2_l:0.2:scop2_r);
xy = [X(:) Y(:)];
p = mvnpdf(xy,ava,sigma);                     %联合概率密度
P = reshape(p,size(X));
 
mesh(X,Y,P);                                  %3-D概率密度图形
name1 = ['二维正态分布 N(',num2str(u1),',',num2str(u2)];
name2 = [',',num2str(sig1),',',num2str(sig2),',',num2str(rho),')'];
name = [name1,name2];
title(name);
end