多维高斯随机变量的仿真

多维高斯随机变量

问题

在进行科学研究的时候，可能会涉及到如何采样多维高斯随机变量的问题，即给定随机变量$\xi \sim \mathcal{N}(0,\mathbf{K}),\xi \in {\mathbb{R}}^n$，得到样本数据${\mathbf{x}}_i,i=1,2,...Q$。这里假设高斯随机变量的均值为零，对于非零的均值，只需通过对样本数据添加一个常量即可。多维高斯随机变量的概率密度函数由均值和协方差矩阵唯一确定。

解决思路

1. 先采样n维标准正态分布得到数据样本${\mathbf{y}}_i$，这可以通过以下方式实现：采样一维标准正态分布，再将每n个一维样本级联成一个n维样本。
2. 对$\mathbf{K}$进行cholesky分解，即$\mathbf{K}=\mathbf{A}{\mathbf{A}}^T$。
3. 获得最终的样本${\mathbf{x}}_i=\mathbf{A}{\mathbf{y}}_i$。

Matlab代码实现

function [ X ] = joint_gaussian( K , num)
% sampling the multi-dimensional Gaussian random variablbe
% K is the covariance matrix
% num is the number of the samples
if K'~=K
   error('Kz must be symmetry');
end
N=size(K,1);
L=chol(K,'lower');
z=randn(N,num);
X=L*z;
end

扩展

循环对称多维复高斯随机变量

在通信仿真中，常常涉及到循环对称多维复高斯随机变量的采样，这里先解释一下什么是循环对称多维复高斯随机变量？它是多维复高斯随机变量的特例。和多维高斯随机变量不同，多维复随机高斯变量的概率密度函数仅由协方差矩阵${\mathbf{K}}_z=E[\xi {\xi }^H]$不能完全确定，还需要伪协方差矩阵${\mathbf{M}}_z=E[\xi {\xi }^T]$共同确定。但如果是循环对称多维复高斯，则其概率密度函数由协方差矩阵完全确定[1]。因为循环对称多维复高斯的伪协方差矩阵为零矩阵。
采样的方法如下所示：

1. 将n维的复随机变量$\xi$转化为2n维的实随机变量$\eta =[Real(\xi ),Imag(\xi )]$。
2. 计算$\eta$的协方差矩阵${\mathbf{K}}_{\eta }=\left[\begin{array}{cc}0.5Real({\mathbf{K}}_z)& -0.5Imag({\mathbf{K}}_z)\\ 0.5Imag({\mathbf{K}}_z)& 0.5Real({\mathbf{K}}_z)\end{array}\right]$
3. 使用上一节的采样方法进行采样，再将其组合成复随机变量。

Matlab代码实现

function [ X ] = joint_complex_gaussian( Kz, num )
% sampling the circulary-symmetric Gaussian random vector
% Kz is the covariance matrix
% num is the number of the samples
if Kz'~=Kz
   error('Kz must be hermitian');
end
Kv=[0.5*real(Kz),-0.5*imag(Kz);0.5*imag(Kz),0.5*real(Kz)];
N=size(Kv,1);
L=chol(Kv,'lower');
z=randn(N,num);
X=L*z;
X=X(1:N/2,:)+1i*X(N/2+1:end,:);
end

代码测试

对采样到的样本数据计算协方差矩阵，再与理论值进行比较，使用相对误差衡量效果，在样本数据为10万时，相对误差在0.5%左右。

clear all;
n = 6;
num = 100000;
% joint gaussian
K = randn(n,n);
K = K*K';
x = joint_gaussian(K,num);
% validate the covariance matrix by relative error
norm(K-x*x'/num)/norm(K)

% joint complex gaussian
Kz = randn(n,n)+1i*randn(n,n);
Kz = Kz*Kz';
x_complex = joint_complex_gaussian(Kz,num);
% validate the covariance matrix by relative error
norm(Kz-x_complex*x_complex'/num)/norm(Kz)

本文使用的符号

${\mathbf{A}}^T$	矩阵转置
${\mathbf{A}}^H$	矩阵共轭转置
$E[]$	数学期望
$Real()$	取复数的实部
$Imag()$	取复数的虚部

参考文献

[1] Gallager R G. Circularly-Symmetric Gaussian random vectors[J]. preprint, 2008:1-9.