在多元正态分布、均匀分布和经验分布中实现拉丁超立方体采样的程序

在多元正态分布、均匀分布和经验分布中实现拉丁超立方体采样的程序。

拉丁超立方体采样概述

拉丁超立方体采样(LHS)是一种分层随机采样方法，能够在较少样本情况下更好地覆盖参数空间。

完整的LHS采样工具箱

1. 核心LHS生成函数

function samples = lhsdesign_custom(n, d, varargin)
% 自定义拉丁超立方体采样
% 输入：
%   n - 样本数量
%   d - 维度
%   varargin - 可选参数：'iterations', 'criterion'
% 输出：
%   samples - n×d的LHS样本矩阵，范围[0,1]^d

    p = inputParser;
    addOptional(p, 'iterations', 5, @isnumeric);
    addOptional(p, 'criterion', 'maximin', @ischar);
    parse(p, varargin{:});
    
    iterations = p.Results.iterations;
    criterion = p.Results.criterion;
    
    best_design = [];
    best_criterion_value = -inf;
    
    for iter = 1:iterations
        % 生成拉丁超立方体设计
        X = lhsdesign_base(n, d);
        
        % 根据准则评估
        switch criterion
            case 'maximin'
                current_value = calculate_maximin_distance(X);
            case 'correlation'
                current_value = -calculate_correlation(X);
            case 'centered'
                current_value = calculate_centered_criterion(X);
            otherwise
                current_value = calculate_maximin_distance(X);
        end
        
        % 选择最佳设计
        if current_value > best_criterion_value
            best_criterion_value = current_value;
            best_design = X;
        end
    end
    
    samples = best_design;
end

function X = lhsdesign_base(n, d)
% 基础LHS生成
    X = zeros(n, d);
    for i = 1:d
        X(:, i) = (randperm(n)' - rand(n, 1)) / n;
    end
end

function d = calculate_maximin_distance(X)
% 计算最大最小距离准则
    n = size(X, 1);
    D = pdist(X);
    d = min(D);
end

function c = calculate_correlation(X)
% 计算相关性准则（绝对值之和）
    R = corr(X);
    c = sum(abs(R(:))) - size(X, 2); % 减去对角线
end

2. 多元正态分布LHS采样

function samples = lhs_norm(n, mu, Sigma, varargin)
% 多元正态分布的拉丁超立方体采样
% 输入：
%   n - 样本数量
%   mu - 均值向量 (d×1)
%   Sigma - 协方差矩阵 (d×d)
%   varargin - LHS参数
% 输出：
%   samples - 来自多元正态分布的LHS样本

    d = length(mu);
    
    % 参数验证
    if size(Sigma, 1) ~= d || size(Sigma, 2) ~= d
        error('协方差矩阵维度必须与均值向量匹配');
    end
    
    % 生成单位LHS样本 [0,1]^d
    lhs_unit = lhsdesign_custom(n, d, varargin{:});
    
    % 转换为标准正态分布
    samples_norm = norminv(lhs_unit, 0, 1);
    
    % 应用Cholesky分解进行相关结构变换
    try
        L = chol(Sigma, 'lower');
    catch
        % 如果协方差矩阵不是正定的，使用最近正定矩阵
        Sigma = nearestSPD(Sigma);
        L = chol(Sigma, 'lower');
    end
    
    % 变换到目标多元正态分布
    samples = bsxfun(@plus, samples_norm * L', mu');
    
    % 验证结果
    fprintf('多元正态分布LHS采样完成:\n');
    fprintf('  样本数: %d, 维度: %d\n', n, d);
    fprintf('  理论均值: %s\n', mat2str(mu', 2));
    fprintf('  样本均值: %s\n', mat2str(mean(samples)', 2));
end

function A = nearestSPD(A)
% 找到最近的对称正定矩阵
    [V, D] = eig((A + A')/2);
    d = diag(D);
    d(d <= 0) = eps;
    A = V * diag(d) * V';
end

3. 多元均匀分布LHS采样

function samples = lhs_uniform(n, lower_bounds, upper_bounds, varargin)
% 多元均匀分布的拉丁超立方体采样
% 输入：
%   n - 样本数量
%   lower_bounds - 下界向量 (1×d)
%   upper_bounds - 上界向量 (1×d)
%   varargin - LHS参数
% 输出：
%   samples - 来自多元均匀分布的LHS样本

    d = length(lower_bounds);
    
    if length(upper_bounds) ~= d
        error('上下界向量维度必须一致');
    end
    
    % 生成单位LHS样本
    lhs_unit = lhsdesign_custom(n, d, varargin{:});
    
    % 线性变换到指定范围
    samples = zeros(n, d);
    for i = 1:d
        samples(:, i) = lower_bounds(i) + ...
                       lhs_unit(:, i) * (upper_bounds(i) - lower_bounds(i));
    end
    
    % 验证结果
    fprintf('多元均匀分布LHS采样完成:\n');
    fprintf('  样本数: %d, 维度: %d\n', n, d);
    fprintf('  参数范围: %s to %s\n', mat2str(lower_bounds, 2), ...
            mat2str(upper_bounds, 2));
end

4. 经验分布LHS采样

function samples = lhs_empirical(n, data, varargin)
% 经验分布的拉丁超立方体采样
% 输入：
%   n - 样本数量
%   data - 经验数据矩阵 (m×d)，每列是一个变量
%   varargin - LHS参数
% 输出：
%   samples - 来自经验分布的LHS样本

    [m, d] = size(data);
    
    if m < n
        warning('样本数量超过数据量，将使用bootstrap补充');
    end
    
    % 生成单位LHS样本
    lhs_unit = lhsdesign_custom(n, d, varargin{:});
    
    % 使用经验逆CDF转换
    samples = zeros(n, d);
    for i = 1:d
        % 对每个维度使用经验分布函数
        sorted_data = sort(data(:, i));
        
        % 经验CDF
        ecdf_vals = (1:m) / m;
        
        % 线性插值得到分位数
        samples(:, i) = interp1(ecdf_vals, sorted_data, lhs_unit(:, i), ...
                               'linear', 'extrap');
    end
    
    % 验证结果
    fprintf('经验分布LHS采样完成:\n');
    fprintf('  样本数: %d, 维度: %d, 原始数据量: %d\n', n, d, m);
    
    % 比较统计量
    for i = 1:d
        orig_mean = mean(data(:, i));
        samp_mean = mean(samples(:, i));
        fprintf('  维度%d - 原始均值: %.4f, 样本均值: %.4f\n', ...
                i, orig_mean, samp_mean);
    end
end

5. 可视化与分析工具

function plot_lhs_comparison(samples, distribution_type, var_names)
% 绘制LHS采样结果的可视化
% 输入：
%   samples - 采样数据
%   distribution_type - 分布类型字符串
%   var_names - 变量名称细胞数组

    [n, d] = size(samples);
    
    if nargin < 3
        var_names = arrayfun(@(i) sprintf('X%d', i), 1:d, 'UniformOutput', false);
    end
    
    figure;
    
    % 散点图矩阵
    if d >= 2
        subplot(2,2,1);
        if d == 2
            scatter(samples(:,1), samples(:,2), 30, 'filled', 'b');
            xlabel(var_names{1}); ylabel(var_names{2});
            grid on;
        else
            plotmatrix(samples);
        end
        title(sprintf('%s LHS采样 - 散点图', distribution_type));
    end
    
    % 边际分布直方图
    subplot(2,2,2);
    if d == 1
        histogram(samples, 20, 'Normalization', 'pdf');
        xlabel(var_names{1}); ylabel('密度');
    else
        histogram(samples(:,1), 20, 'Normalization', 'pdf');
        xlabel(var_names{1}); ylabel('密度');
    end
    title('边际分布');
    
    % 空间填充性评估
    subplot(2,2,3);
    if d >= 2
        % 计算最小距离
        D = pdist(samples);
        min_distances = zeros(n,1);
        for i = 1:n
            distances = sqrt(sum((samples - samples(i,:)).^2, 2));
            distances(i) = inf; % 排除自身
            min_distances(i) = min(distances);
        end
        
        plot(1:n, sort(min_distances), 'b-o', 'LineWidth', 1.5);
        xlabel('样本索引'); ylabel('最小距离');
        title('空间填充性分析');
        grid on;
    end
    
    % 相关性矩阵
    subplot(2,2,4);
    if d >= 2
        R = corr(samples);
        imagesc(R);
        colorbar;
        title('样本相关性矩阵');
        set(gca, 'XTick', 1:d, 'XTickLabel', var_names);
        set(gca, 'YTick', 1:d, 'YTickLabel', var_names);
    end
    
    sgtitle(sprintf('%s分布LHS采样分析 (n=%d, d=%d)', ...
            distribution_type, n, d));
end

function [stats] = analyze_lhs_performance(samples, true_params)
% 分析LHS采样性能
% 输入：
%   samples - LHS样本
%   true_params - 真实参数（可选）
    
    stats = struct();
    [n, d] = size(samples);
    
    % 基本统计量
    stats.sample_size = n;
    stats.dimensions = d;
    stats.sample_mean = mean(samples)';
    stats.sample_std = std(samples)';
    stats.sample_corr = corr(samples);
    
    % 空间填充性指标
    D = pdist(samples);
    stats.min_distance = min(D);
    stats.avg_distance = mean(D);
    stats.space_filling = stats.min_distance;
    
    % 与真实参数比较（如果提供）
    if nargin > 1
        if isfield(true_params, 'mean')
            stats.mean_error = norm(stats.sample_mean - true_params.mean);
        end
        if isfield(true_params, 'cov')
            stats.cov_error = norm(cov(samples) - true_params.cov, 'fro');
        end
    end
    
    % 显示结果
    fprintf('\n=== LHS采样性能分析 ===\n');
    fprintf('样本数: %d, 维度: %d\n', n, d);
    fprintf('最小样本距离: %.4f\n', stats.min_distance);
    fprintf('平均样本距离: %.4f\n', stats.avg_distance);
    
    if isfield(stats, 'mean_error')
        fprintf('均值估计误差: %.4f\n', stats.mean_error);
    end
end

使用

示例1：多元正态分布LHS采样

% 定义二元正态分布参数
mu = [2, 5];                    % 均值
Sigma = [1, 0.6; 0.6, 2];       % 协方差矩阵

% 生成LHS样本
n_samples = 100;
samples_norm = lhs_norm(n_samples, mu, Sigma, 'iterations', 10, 'criterion', 'maximin');

% 可视化与分析
plot_lhs_comparison(samples_norm, '多元正态', {'X1', 'X2'});
stats_norm = analyze_lhs_performance(samples_norm, struct('mean', mu', 'cov', Sigma));

示例2：多元均匀分布LHS采样

% 定义均匀分布参数
lower_bounds = [0, -1, 2];      % 下界
upper_bounds = [5, 3, 8];       % 上界

% 生成LHS样本
samples_unif = lhs_uniform(50, lower_bounds, upper_bounds, 'criterion', 'centered');

% 可视化
plot_lhs_comparison(samples_unif, '多元均匀', {'U1', 'U2', 'U3'});

示例3：经验分布LHS采样

% 生成示例经验数据（金融收益）
rng(42);
historical_returns = mvnrnd([0.08, 0.12], [0.2, 0.1; 0.1, 0.3], 1000);

% 基于历史数据的LHS采样
samples_empirical = lhs_empirical(200, historical_returns);

% 比较原始数据与LHS样本
figure;
subplot(1,2,1);
scatter(historical_returns(:,1), historical_returns(:,2), 20, 'b', 'filled');
title('原始历史数据'); xlabel('资产1收益'); ylabel('资产2收益');

subplot(1,2,2);
scatter(samples_empirical(:,1), samples_empirical(:,2), 40, 'r', 'filled');
title('LHS采样数据'); xlabel('资产1收益'); ylabel('资产2收益');

示例4：综合比较

% 比较不同采样方法
n = 50;
d = 2;

% 随机采样
random_samples = rand(n, d);

% LHS采样
lhs_samples = lhsdesign_custom(n, d, 'iterations', 5);

% 可视化比较
figure;
subplot(1,2,1);
scatter(random_samples(:,1), random_samples(:,2), 40, 'b', 'filled');
title('随机采样'); xlabel('X1'); ylabel('X2'); grid on;

subplot(1,2,2);
scatter(lhs_samples(:,1), lhs_samples(:,2), 40, 'r', 'filled');
title('LHS采样'); xlabel('X1'); ylabel('X2'); grid on;

% 计算性能指标
fprintf('随机采样最小距离: %.4f\n', min(pdist(random_samples)));
fprintf('LHS采样最小距离: %.4f\n', min(pdist(lhs_samples)));

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

参数调优建议

1. 样本数量：通常为维度数的10-50倍

2. 迭代次数：5-20次以获得较好的空间填充性

3. 选择准则：

   • maximin：最大化最小距离（推荐）
   • correlation：最小化相关性
   • centered：样本更靠近区间中心