【发布时间】:2017-04-18 19:22:03
【问题描述】:
我想知道如何在 matlab 中绘制样本,其中我有精度矩阵和均值作为输入参数。
我知道 mvnrnd 是一种典型的方法,但它需要协方差矩阵(即精度的倒数)作为参数。
我只有精度矩阵,由于计算问题,我无法反转我的精度矩阵,因为它需要太长时间(我的尺寸约为 2000*2000)
【问题讨论】:
标签: matlab matrix probability
我想知道如何在 matlab 中绘制样本,其中我有精度矩阵和均值作为输入参数。
我知道 mvnrnd 是一种典型的方法,但它需要协方差矩阵(即精度的倒数)作为参数。
我只有精度矩阵,由于计算问题,我无法反转我的精度矩阵,因为它需要太长时间(我的尺寸约为 2000*2000)
【问题讨论】:
标签: matlab matrix probability
好问题。请注意,您可以通过the relevant Wikipedia article 中描述的过程使用来自标准正态分布的样本从多元正态分布生成样本。
基本上,这归结为评估A*z + mu,其中z 是从标准正态分布中采样的独立随机变量向量,mu 是均值向量,A*A' = Sigma 是协方差矩阵。由于您有后一个数量的倒数,即inv(Sigma),您可能可以进行 Cholesky 分解(参见chol)来确定A 的倒数。然后您需要评估A * z。如果您只知道inv(A),这仍然可以在不执行矩阵求逆的情况下通过求解线性系统(例如通过反斜杠运算符)来完成。
Cholesky 分解对您来说可能仍然存在问题,但我希望这会有所帮助。
【讨论】:
如果你想从 N(μ,Q-1) 中采样并且只有 Q 可用,你可以对 Q, L 进行 Cholesky 分解,使得 LLT支持>=Q。接下来取 LT、L-T 的逆,并从标准正态分布 N(0, I) 中采样 Z。
考虑到L-T是上三角dxd矩阵,Z是d维列向量, μ + L-TZ 将分布为 N(μ, Q-1)。
如果您希望避免取 L 的倒数,您可以改为通过反向代入求解三角方程组 LTv=Z。 μ+v 将分布为 N(μ, Q-1)。
一些说明性的matlab代码:
% 制作一个 2x2 协方差矩阵和一个均值向量
covm = [3 0.4*(sqrt(3*7)); 0.4*(sqrt(3*7)) 7];
亩 = [100; 2];
% 获取精度矩阵
Q = inv(covm);
%对Q进行Cholesky分解(matlab中的chol已经返回上三角因子)
L = chol(Q);
%从标准二元正态分布中抽取 2000 个样本
Z = normrnd(0,1, [2, 2000]);
%求解系统并添加均值
X = repmat(mu, 1, 2000)+L\Z;
%检查结果
平均(X')
var(X')
corrcoef(X')
% 与协方差矩阵中的采样进行比较
Y=mvnrnd(mu,covm, 2000)';
平均(Y')
var(Y')
corrcoef(Y')
scatter(X(1,:), X(2,:),'b')
等一下
scatter(Y(1,:), Y(2,:), 'r')
为了提高效率,我想你可以搜索一些可以有效解决三角系统的包。
【讨论】: