%% Sobol敏感性分析实现 - 修正版本
% 处理58行10列数据集,前7列为输入变量,后3列为输出变量
% 使用正确的Sobol敏感性分析方法
clear; clc; close all;
%% 主函数
fprintf('Sobol敏感性分析程序启动...\n');
fprintf('==================================================\n');
try
% 1. 加载数据
data = load_user_data();
fprintf('\n数据加载完成: %d行 × %d列\n', size(data,1), size(data,2));
fprintf('数据预览:\n');
% 创建变量名
varNames = [strseq('X',1:7), strseq('Y',1:3)];
dataTable = array2table(data, 'VariableNames', varNames);
disp(head(dataTable));
% 2. 准备分析数据
[X, Y, problem] = prepare_sobol_analysis(data);
% 3. 执行真正的Sobol分析
results = perform_true_sobol_analysis(X, Y, problem);
% 4. 识别高敏感性变量
high_sensitivity_vars = identify_high_sensitivity_variables(results, 0.1);
% 5. 可视化结果
visualize_sobol_results(results);
fprintf('\n==================================================\n');
fprintf('分析完成!\n');
fprintf('==================================================\n');
catch ME
fprintf('分析过程中出现错误: %s\n', ME.message);
rethrow(ME);
end
%% 辅助函数:生成序列名称
function names = strseq(prefix, indices)
names = arrayfun(@(i) sprintf('%s%d', prefix, i), indices, 'UniformOutput', false);
end
%% 加载用户数据
function data = load_user_data()
fprintf('正在加载数据集...\n');
try
data = readmatrix('分析数据.csv');
fprintf('成功从 筛选后的数据.xlsx 加载数据\n');
catch
fprintf('无法加载数据文件,请检查文件是否存在或格式是否正确\n');
fprintf('将使用示例数据进行演示\n');
data = generate_sample_data();
end
% 验证数据维度
if size(data, 2) ~= 10
fprintf('警告: 数据列数不是10,当前为%d列\n', size(data, 2));
fprintf('请确保数据格式为:前7列输入,后3列输出\n');
end
% 检查数据中是否有NaN或Inf值
if any(isnan(data(:)) | isinf(data(:)))
fprintf('警告: 数据中包含NaN或Inf值,将尝试处理...\n');
% 用列均值替换NaN
for i = 1:size(data, 2)
col_data = data(:, i);
nan_indices = isnan(col_data);
if any(nan_indices)
col_mean = mean(col_data(~nan_indices & ~isinf(col_data)));
data(nan_indices, i) = col_mean;
fprintf(' 列 %d: 替换了 %d 个NaN值\n', i, sum(nan_indices));
end
% 处理Inf值
inf_indices = isinf(col_data);
if any(inf_indices)
col_median = median(col_data(~inf_indices & ~isnan(col_data)));
data(inf_indices, i) = col_median;
fprintf(' 列 %d: 替换了 %d 个Inf值\n', i, sum(inf_indices));
end
end
end
end
%% 生成示例数据(仅在无法加载用户数据时使用)
function data = generate_sample_data()
fprintf('使用示例数据进行演示...\n');
% 生成示例数据集:58行10列,前7列输入,后3列输出
rng(42); % 设置随机种子
% 生成输入变量 (X1-X7)
n_samples = 58;
n_inputs = 7;
% 创建输入变量 - 使用均匀分布
X = rand(n_samples, n_inputs);
% 生成输出变量 (Y1-Y3),基于输入变量的非线性组合,包含交互效应
Y1 = 2.0 * X(:, 1) + 1.5 * X(:, 2).^2 + 0.5 * X(:, 3) .* X(:, 4) + 0.1 * randn(n_samples, 1);
Y2 = 1.0 * X(:, 2) + 3.0 * X(:, 5).^3 + 0.8 * X(:, 6) + 0.2 * X(:, 1) .* X(:, 7) + 0.1 * randn(n_samples, 1);
Y3 = 0.5 * X(:, 1) .* X(:, 7) + 2.0 * X(:, 2) + 1.0 * X(:, 6).^2 + 0.3 * X(:, 3) .* X(:, 5) + 0.1 * randn(n_samples, 1);
% 组合数据
data = [X, Y1, Y2, Y3];
end
%% 准备Sobol分析
function [X, Y, problem] = prepare_sobol_analysis(data)
% 提取输入和输出
X = data(:, 1:7); % 前7列是输入
Y = data(:, 8:10); % 后3列是输出
% 检查输出变量是否有零方差
for i = 1:size(Y, 2)
if var(Y(:, i)) == 0
fprintf('警告: 输出变量 Y%d 方差为零,添加微小噪声...\n', i);
Y(:, i) = Y(:, i) + 1e-10 * randn(size(Y, 1), 1);
end
end
% 计算每个输入变量的实际范围
bounds = [min(X); max(X)]';
% 定义问题
problem = struct();
problem.num_vars = 7;
problem.names = strseq('X', 1:7);
problem.bounds = bounds; % 使用实际数据的范围
fprintf('输入变量范围:\n');
for i = 1:7
fprintf(' %s: [%.4f, %.4f]\n', problem.names{i}, bounds(i,1), bounds(i,2));
end
end
%% 执行真正的Sobol分析
function results = perform_true_sobol_analysis(X, Y, problem)
fprintf('\n正在执行真正的Sobol敏感性分析...\n');
results = struct();
% 对每个输出变量进行分析
for i = 1:size(Y, 2)
output_name = sprintf('Y%d', i);
fprintf('\n分析输出变量: %s\n', output_name);
y = Y(:, i);
% 使用基于方差的Sobol方法
try
% 方法1: 使用Saltelli方法近似计算Sobol指数
[S1, ST] = calculate_sobol_indices(X, y);
catch ME
fprintf('Sobol分析失败: %s\n', ME.message);
fprintf('使用替代方法...\n');
% 方法2: 使用基于回归的方法
[S1, ST] = calculate_regression_based_indices(X, y);
end
% 创建结果结构
results.(output_name) = struct();
results.(output_name).S1 = S1; % 一阶敏感性指数
results.(output_name).ST = ST; % 总效应指数
results.(output_name).names = problem.names;
% 打印结果
fprintf('\n%s的Sobol敏感性分析结果:\n', output_name);
fprintf('变量\t\t一阶指数\t总效应指数\t差值(交互效应)\n');
for j = 1:length(problem.names)
interaction_effect = ST(j) - S1(j);
fprintf('%s\t\t%.4f\t\t%.4f\t\t%.4f\n', problem.names{j}, S1(j), ST(j), interaction_effect);
end
% 计算并显示交互效应占比
total_interaction = sum(ST - S1);
fprintf('总交互效应: %.4f\n', total_interaction);
end
end
%% 计算Sobol指数 - 基于Saltelli方法
function [S1, ST] = calculate_sobol_indices(X, y)
n_vars = size(X, 2);
n_samples = size(X, 1);
% 初始化结果
S1 = zeros(n_vars, 1);
ST = zeros(n_vars, 1);
% 计算总方差
total_variance = var(y);
if total_variance == 0
error('输出变量方差为零,无法计算Sobol指数');
end
% 对每个变量计算一阶效应和总效应
for i = 1:n_vars
% 一阶效应:固定其他变量,只变化第i个变量
% 使用条件方差估计
% 分组数据:按第i个变量的值分组
[sorted_X, sort_idx] = sort(X(:, i));
sorted_y = y(sort_idx);
% 使用移动窗口计算条件期望
window_size = max(3, floor(n_samples / 10)); % 窗口大小
conditional_means = zeros(n_samples, 1);
for j = 1:n_samples
start_idx = max(1, j - floor(window_size/2));
end_idx = min(n_samples, j + floor(window_size/2));
conditional_means(j) = mean(sorted_y(start_idx:end_idx));
end
% 计算一阶效应方差
variance_S1 = var(conditional_means);
S1(i) = variance_S1 / total_variance;
% 总效应:使用基于抽样的方法
% 固定第i个变量,随机重排其他变量
n_permutations = min(100, n_samples);
total_effect_variances = zeros(n_permutations, 1);
for p = 1:n_permutations
% 创建新样本:固定第i个变量,重排其他变量
X_permuted = X;
other_vars = setdiff(1:n_vars, i);
% 重排其他变量
perm_idx = randperm(n_samples);
X_permuted(:, other_vars) = X(perm_idx, other_vars);
% 使用KNN预测输出
if n_samples > 10
k = min(5, floor(n_samples/10));
y_pred = knn_predict(X, y, X_permuted, k);
total_effect_variances(p) = var(y_pred);
else
% 样本太少,使用简单方法
total_effect_variances(p) = total_variance;
end
end
% 计算总效应
ST(i) = mean(total_effect_variances) / total_variance;
end
% 确保指数在合理范围内
S1 = max(0, min(1, S1));
ST = max(S1, min(1, ST)); % 总效应应该不小于一阶效应
end
%% 使用KNN进行预测
function y_pred = knn_predict(X_train, y_train, X_test, k)
n_test = size(X_test, 1);
y_pred = zeros(n_test, 1);
for i = 1:n_test
% 计算距离
distances = sqrt(sum((X_train - X_test(i, :)).^2, 2));
% 找到最近的k个邻居
[~, idx] = sort(distances);
nearest_indices = idx(1:min(k, length(idx)));
% 使用邻居的平均值作为预测
y_pred(i) = mean(y_train(nearest_indices));
end
end
%% 基于回归的方法计算敏感性指数
function [S1, ST] = calculate_regression_based_indices(X, y)
n_vars = size(X, 2);
% 标准化数据
X_scaled = zscore(X);
y_scaled = zscore(y);
% 拟合线性模型获取一阶效应
mdl = fitlm(X_scaled, y_scaled);
coefficients = mdl.Coefficients.Estimate(2:end); % 排除截距
% 一阶效应:标准化回归系数的平方
S1 = (coefficients .^ 2) / sum(coefficients .^ 2);
% 总效应:使用随机森林重要性作为近似
try
% 使用随机森林估计总效应
tree = fitrensemble(X_scaled, y_scaled, 'Method', 'Bag', 'NumLearningCycles', 50);
imp = oobPermutedPredictorImportance(tree);
ST = imp / sum(imp);
catch
% 如果随机森林失败,使用线性模型加上交互项
ST = calculate_total_effect_with_interactions(X_scaled, y_scaled);
end
% 确保总效应不小于一阶效应
ST = max(S1, ST);
end
%% 计算包含交互效应的总效应
function ST = calculate_total_effect_with_interactions(X, y)
n_vars = size(X, 2);
ST = zeros(n_vars, 1);
for i = 1:n_vars
% 创建包含所有变量和与第i个变量交互项的模型
interaction_terms = [];
for j = 1:n_vars
if j ~= i
interaction_terms = [interaction_terms, X(:, i) .* X(:, j)];
end
end
% 组合特征
features = [X, interaction_terms];
% 拟合模型
mdl = fitlm(features, y);
% 计算第i个变量及其交互项的贡献
var_contrib = 0;
for j = 1:length(mdl.CoefficientNames)
coef_name = mdl.CoefficientNames{j};
if contains(coef_name, sprintf('x%d', i)) || ...
(contains(coef_name, 'x') && contains(coef_name, sprintf('x%d_', i)))
var_contrib = var_contrib + mdl.Coefficients.Estimate(j)^2;
end
end
ST(i) = var_contrib;
end
% 归一化
ST = ST / sum(ST);
end
%% 识别高敏感性变量
function high_sensitivity_vars = identify_high_sensitivity_variables(results, threshold)
if nargin < 2
threshold = 0.1;
end
fprintf('\n%s\n', repmat('=', 1, 60));
fprintf('高敏感性变量识别结果:\n');
fprintf('%s\n', repmat('=', 1, 60));
high_sensitivity_vars = struct();
output_names = fieldnames(results);
for i = 1:length(output_names)
output = output_names{i};
result = results.(output);
fprintf('\n%s的高敏感性变量 (阈值: %.1f):\n', output, threshold);
% 基于总效应指数识别高敏感性变量
high_indices = find(result.ST > threshold & ~isnan(result.ST));
if ~isempty(high_indices)
high_sensitivity_vars.(output) = {};
for j = 1:length(high_indices)
idx = high_indices(j);
var_name = result.names{idx};
S1 = result.S1(idx);
ST = result.ST(idx);
interaction = ST - S1;
fprintf(' %s: 一阶=%.4f, 总效应=%.4f, 交互效应=%.4f\n', ...
var_name, S1, ST, interaction);
high_sensitivity_vars.(output){end+1} = {var_name, S1, ST, interaction};
end
else
fprintf(' 无变量超过阈值\n');
end
end
end
%% 可视化Sobol结果
function visualize_sobol_results(results)
fprintf('\n正在生成Sobol敏感性分析图表...\n');
output_names = fieldnames(results);
n_outputs = length(output_names);
% 创建主图
fig = figure('Position', [100, 100, 1400, 1000]);
for i = 1:n_outputs
output = output_names{i};
result = results.(output);
subplot(2, 2, i);
% 准备数据
variables = result.names;
s1_values = result.S1;
st_values = result.ST;
interaction_values = st_values - s1_values;
% 检查是否有NaN值
valid_indices = ~isnan(s1_values) & ~isnan(st_values);
if ~all(valid_indices)
fprintf('警告: %s 中有NaN值,将跳过这些变量\n', output);
variables = variables(valid_indices);
s1_values = s1_values(valid_indices);
st_values = st_values(valid_indices);
interaction_values = interaction_values(valid_indices);
end
if isempty(variables)
text(0.5, 0.5, '无有效数据', 'HorizontalAlignment', 'center', 'FontSize', 14);
title(sprintf('%s - 无有效数据', output), 'FontSize', 12);
continue;
end
% 创建堆叠柱状图显示一阶效应和交互效应
x = 1:length(variables);
% 绘制一阶效应
bars1 = bar(x, s1_values, 0.8, 'FaceColor', [0.2, 0.6, 0.8]);
hold on;
% 绘制交互效应(堆叠在一阶效应之上)
bars2 = bar(x, interaction_values, 0.8, 'FaceColor', [0.8, 0.4, 0.2]);
% 设置堆叠
bars2.BaseValue = 0;
% 添加总效应线
plot(x, st_values, 'ko-', 'LineWidth', 2, 'MarkerSize', 6, 'MarkerFaceColor', 'k');
xlabel('输入变量');
ylabel('敏感性指数');
title(sprintf('%s的Sobol敏感性分析', output), 'FontSize', 12);
set(gca, 'XTick', x, 'XTickLabel', variables);
xtickangle(45);
% legend({'一阶效应 (S1)', '交互效应', '总效应 (ST)'}, 'Location', 'best');
grid on;
% 添加阈值线
yline(0.1, '--r', '阈值=0.1', 'LabelHorizontalAlignment', 'left', 'FontSize', 8);
% 添加数值标签
for j = 1:length(x)
text(x(j), st_values(j) + 0.02, sprintf('ST=%.3f', st_values(j)), ...
'HorizontalAlignment', 'center', 'VerticalAlignment', 'bottom', 'FontSize', 8);
end
end
% 添加综合比较图在第四个位置
if n_outputs >= 3
subplot(2, 2, 4);
% 准备比较数据
s1_matrix = zeros(n_outputs, length(results.(output_names{1}).names));
st_matrix = zeros(n_outputs, length(results.(output_names{1}).names));
for i = 1:n_outputs
s1_matrix(i, :) = results.(output_names{i}).S1;
st_matrix(i, :) = results.(output_names{i}).ST;
end
% 检查是否有NaN值
if any(isnan(s1_matrix(:))) || any(isnan(st_matrix(:)))
fprintf('警告: 比较图数据中包含NaN值,将用0替换\n');
s1_matrix(isnan(s1_matrix)) = 0;
st_matrix(isnan(st_matrix)) = 0;
end
% 创建分组柱状图
x = 1:size(s1_matrix, 2);
width = 0.35;
for i = 1:n_outputs
bar_positions = x + (i-1) * width - (n_outputs-1)*width/2;
bar(bar_positions, st_matrix(i, :), width, 'FaceAlpha', 0.7);
hold on;
end
xlabel('输入变量');
ylabel('总敏感度指数 (ST)');
title('各输出变量的输入变量总敏感性比较', 'FontSize', 12);
set(gca, 'XTick', x, 'XTickLabel', results.(output_names{1}).names);
% legend(output_names, 'Location', 'best');
grid on;
end
% 调整布局
sgtitle('Sobol敏感性分析结果 (S1: 一阶效应, ST: 总效应)', 'FontSize', 16, 'FontWeight', 'bold');
% 保存图表
saveas(fig, 'sobol_sensitivity_analysis_corrected.png');
fprintf('图表已保存为: sobol_sensitivity_analysis_corrected.png\n');
end上面的代码运行之后得到的结果我认为存在问题,得到的一阶灵敏度指标与全局灵敏度指标相差值过大,且得到的结果准确性我认为不强,请帮我对代码进行修改完善,保证运行之后得到的结果更准确
最新发布