Matlab拉丁超立方抽样（与Cholesky 分解法重排序）

魔法帽

已于 2024-03-25 16:37:33 修改

阅读量1.9k

点赞数 22

文章标签：概率论 matlab 算法

于 2024-03-25 16:30:32 首次发布

本文链接：https://blog.csdn.net/Liao_223/article/details/137016908

版权

本文介绍了如何使用Matlab中的lhsnorm和lhsdesign函数对6个随机变量进行抽样，特别关注了正态分布和均匀分布的处理。通过Cholesky分解法降低不同变量列之间的相关性，发现lhsnorm在输出时已进行了类似处理，证实了无需额外重排。

摘要由CSDN通过智能技术生成

简介：

对6个随机变量进行拉丁超立方抽样：

3个正态分布随机变量： ${{x}_{1}}\sim N(0,{{0.05}^{2}})$ , ${{x}_{2}}\sim N(0,{{0.05}^{2}})$ , ${{x}_{3}}\sim N(1,{{0.02}^{2}})$

3个均匀分布随机变量： ${{x}_{4},{x}_{5},{x}_{6}}\sim U(0,1)$

超立方抽样代码：

要点：调用matlab中的

lhsnorm()函数进行正态分布抽样，

lhsdesign()函数进行均匀分布抽样。

代码中Xrms为抽样矩阵的相关系数矩阵方均根，用于衡量不同变量列之间的相关性，相关性越低越好。

注意：输出结果已经被重排，相关性较低。

clc;clear;close all;

%% 测试三个正态分布
% 产生抽样矩阵XN（N行K列）
N = 500;                   % 数据条数
mu = [0, 0, 1];            % 正态分布均值
K = length(mu);            % 变量个数
sigma0 = [0.05,0.05,0.02]; % 标准差
sigma = diag(sigma0.^2);   % 协方差矩阵（令协方差值都为0）
XN = lhsnorm(mu, sigma, N);

PQ = corrcoef(XN);         % 相关系数矩阵（KxK）
Xrms = sqrt(sum(PQ.^2,"all")-K)/sqrt(K*(K-1));  % 相关系数矩阵方均根

figure(1)
plot(XN(:,1),XN(:,2),'*');
title('正态分布的LHS采样')

%% 三个均匀分布
XB = lhsdesign(N,3);       % 3列0-1的均匀分布

figure(2)
plot(XB(:,1),XN(:,2),'*');
title('正态分布与均匀分布LHS采样')

输出画图结果

Cholesky 分解法重排序

为了使得抽样结果更加具有代表性，需要降低不同变量列之间的相关性，即对列中的元素进行重排列。方法有：Gram-Schmidt序列正交化方法、Cholesky分解法和随机排列法。

其原理请查看论文：《基于拉丁超立方采样的含风电电力系统的概率可靠性评估》

代码如下：

% 利用Cholesky分解法降低XN的列相关性
% 获得初始排序矩阵QN
QN = zeros(N,K);
for i=1:K
    QN(:,i)=randperm(N);
end

PQ = corrcoef(QN);         % 相关系数矩阵（KxK）
Qrms = sqrt(sum(PQ.^2,"all")-K)/sqrt(K*(K-1));  % 相关系数矩阵方均根
while Qrms>1e-10
    L = chol(PQ,'lower');      % Cholesky分解法得下三角矩阵
    QN = QN/L;                 % 更新排序矩阵
    PQ = corrcoef(QN);         % 相关系数矩阵（KxK）
    Qrms = sqrt(sum(PQ.^2,"all")-K)/sqrt(K*(K-1));  % 相关系数矩阵方均根
    XN = re_arrange(XN,QN);
end

PQ = corrcoef(XN);         % 相关系数矩阵（KxK）
Xrms2 = sqrt(sum(PQ.^2,"all")-K)/sqrt(K*(K-1));  % 相关系数矩阵方均根

function M_AA = re_arrange(MA,MB)
% 按照矩阵MB的大小顺序每一列按行排序MA
    [~,b2]=sort(MB,1,'descend');
    M_AA = zeros(size(MB));
    new_indices = zeros(size(MB));
    for j = 1:size(MB,2)
        for i = 1:size(MB,1)
            new_indices(b2(i,j),j) = i;          % 得到排序序号
        end
    end
    for j = 1:size(MB,2)
        for i = 1:size(MB,1)
            M_AA(new_indices(i,j),j) = MA(i,j);  % 重排
        end
    end
end

结果显示：Xrms在重排序前后变化较小，甚至有所变大，均在0~0.06间浮动。

得出结论：lhsnorm()在输出结果时已经对结果进行类似Cholesky分解法的重排，故无需再进行重排。

验证结论：在 XN = lhsnorm(mu, sigma, N); 代码后面加入XN = sort(XN,1); 此时Xrms=0.9997，进行Cholesky 分解法重排序后，Xrms2=0.0480，证明该Cholesky 分解法有效。