MATLAB实现多分布拉丁超立方采样实用程序

MATLAB实现多分布拉丁超立方采样实用程序，支持多元正态分布、均匀分布和经验分布，并包含相关性控制功能：

---

一、核心代码

function [samples, meta] = multiLHS(nSamples, nVars, distType, varargin)
    % 多分布拉丁超立方采样
    % 输入参数：
    % nSamples - 样本数量
    % nVars - 变量维度
    % distType - 分布类型 ('normal', 'uniform', 'empirical')
    % varargin - 可选参数：
    %   'mu' - 均值向量 (仅normal)
    %   'sigma' - 协方差矩阵 (仅normal)
    %   'bounds' - 边界矩阵 [下限;上限] (仅uniform)
    %   'data' - 经验分布数据 (仅empirical)
    %   'correlation' - 相关性矩阵 (可选)
    % 输出：
    % samples - 采样结果 (nSamples×nVars)
    % meta - 元数据结构体 (包含分布参数)

    % 参数解析
    p = inputParser;
    addParameter(p, 'mu', [], @(x) isnumeric(x) && isvector(x));
    addParameter(p, 'sigma', [], @(x) isnumeric(x) && ismatrix(x));
    addParameter(p, 'bounds', [], @(x) isnumeric(x) && size(x,2)==2);
    addParameter(p, 'data', [], @(x) isnumeric(x));
    addParameter(p, 'correlation', eye(nVars), @(x) isnumeric(x) && ismatrix(x));
    parse(p, varargin{:});
    
    % 基础LHS采样
    baseSamples = lhsdesign(nSamples, nVars);
    
    % 分布转换
    switch lower(distType)
        case 'normal'
            validateattributes(p.Results.mu, {'numeric'}, {'vector', 'numel',nVars});
            validateattributes(p.Results.sigma, {'numeric'}, {'matrix', 'size',[nVars,nVars]});
            samples = mvnrnd(p.Results.mu, p.Results.sigma, nSamples);
            
        case 'uniform'
            validateattributes(p.Results.bounds, {'numeric'}, {'size',[2,nVars]});
            samples = zeros(nSamples,nVars);
            for i = 1:nVars
                samples(:,i) = p.Results.bounds(1,i) + ...
                             (p.Results.bounds(2,i)-p.Results.bounds(1,i)) * baseSamples(:,i);
            end
            
        case 'empirical'
            validateattributes(p.Results.data, {'numeric'}, {'2d'});
            % 核密度估计
            kde = fitgmdist(p.Results.data, 1);
            samples = random(kde, nSamples, nVars);
            
        otherwise
            error('不支持的分布类型');
    end
    
    % 相关性调整
    if ~isequal(correlation, eye(nVars))
        [U,S,V] = svd(correlation);
        samples = samples * U * diag(1./sqrt(diag(S))) * U';
    end
    
    % 元数据记录
    meta = struct(...
        'distribution', distType,...
        'mu', p.Results.mu,...
        'sigma', p.Results.sigma,...
        'bounds', p.Results.bounds,...
        'correlation', correlation...
    );
end

---

二、功能特性

1. 多元正态分布支持

% 示例：生成100个样本，3维正态分布
mu = [1,2,3];
sigma = [2,0.5,0; 0.5,1,0.2; 0,0.2,1];
samples = multiLHS(100,3,'normal','mu',mu,'sigma',sigma);

2. 均匀分布支持

% 示例：生成50个样本，2维均匀分布 [0,1]×[5,10]
bounds = [0,5; 1,10];
samples = multiLHS(50,2,'uniform','bounds',bounds);

3. 经验分布支持

% 示例：基于历史数据生成样本
data = randn(1000,2); % 假设历史数据
samples = multiLHS(200,2,'empirical','data',data);

4. 相关性控制

% 示例：添加相关性矩阵
rho = [1 0.8; 0.8 1];
samples = multiLHS(100,2,'normal','mu',[0,0],'sigma',eye(2),'correlation',rho);

---

三、算法验证

1. 正态分布验证

% 生成样本
mu = [0,0];
sigma = [1,0.5;0.5,1];
samples = multiLHS(1e4,2,'normal','mu',mu,'sigma',sigma);

% 统计验证
mean_val = mean(samples);
cov_val = cov(samples);
disp(['均值误差: ', num2str(norm(mean_val-mu))]);
disp(['协方差误差: ', num2str(norm(cov_val-sigma))]);

2. 相关性验证

% 生成相关样本
rho = [1 0.7; 0.7 1];
samples = multiLHS(1e4,2,'normal','correlation',rho);

% 计算实际相关系数
actual_rho = corr(samples(:,1),samples(:,2));
disp(['实际相关系数: ', num2str(actual_rho)]);

---

四、性能优化

1. 大规模数据处理

% 分块生成 (适用于nSamples>1e5)
chunkSize = 1e4;
nChunks = ceil(nSamples/chunkSize);
samples = zeros(nSamples,nVars);
for i = 1:nChunks
    startIdx = (i-1)*chunkSize +1;
    endIdx = min(i*chunkSize,nSamples);
    samples(startIdx:endIdx,:) = multiLHS(endIdx-startIdx+1,nVars,distType,varargin{:});
end

2. GPU加速

% GPU加速版本
if canUseGPU
    baseSamples = gpuArray(lhsdesign(nSamples,nVars));
    % 后续计算使用gpuArray...
end

参考代码 多元正态分布、均匀分布和经验分布中实现拉丁超立方体采样的采样实用程序 www.youwenfan.com/contentcnm/64716.html

五、应用场景示例

1. 不确定性量化

% 定义参数范围
bounds = [0,10; 1,5; 0.1,2];
% 生成样本
samples = multiLHS(1000,3,'uniform','bounds',bounds);
% 运行仿真模型
results = runModel(samples);
% 计算概率分布
histogram(results(:,1),100);

2. 敏感性分析

% 生成相关样本
rho = [1 0.3; 0.3 1];
samples = multiLHS(500,2,'normal','correlation',rho);
% 计算Sobol指数
S = sobolset(2,'NumSamples',500);
Y = simModel(samples);
S = analyze(SobolIndices,S,Y);

---

六、错误处理机制

% 输入验证示例
try
    samples = multiLHS(-10,3,'normal'); % 错误样本数
catch ME
    disp(['输入错误: ', ME.message]);
end

try
    samples = multiLHS(100,3,'normal','mu',[1,2]); % 维度不匹配
catch ME
    disp(['参数错误: ', ME.message]);
end

---

该实用程序通过模块化设计实现了多分布类型的拉丁超立方采样，支持相关性控制和性能优化，适用于工程仿真、风险评估等场景。用户可根据需要扩展功能模块，如添加自定义分布类型或集成代理模型加速。