第十四届APMCM亚太地区大学生数学建模竞赛B题解题思路深度解析

首先说明思路并不是最优的,但是希望对大家有所启发

第一小问:分析数据并可视化关键指标

第一小问因为问题是询问多个指标中,哪些指标与洪水发生密切相关,哪些指标是不相关,根据表中数据其实很好分析出我们可以直接作图画出用每个参数与洪水发生概率的关系,可以使用用皮尔逊相关系数或者斯皮尔曼等级相关系数来分别找出服从以及不服从的数据。

方法1:使用皮尔逊相关系数(如果数据服从正态分布)

计算每个特征与洪水概率之间的皮尔逊相关系数,将关系可视化便于识别每个系数之间特征。

% 加载数据  
data_train = readtable('train.xlsx');  
  
% 提取特征和标签  
features = data_train{:, 2:end-1};  % 假设最后一列是洪水概率  
labels = data_train.洪水概率;  
  
% 计算皮尔逊相关系数矩阵  
corr_matrix_pearson = corr(features, 'Type', 'Pearson', 'Rows', 'complete');  
  
% 可视化皮尔逊相关系数矩阵  
figure;  
imagesc(corr_matrix_pearson);  
colorbar;  
xlabel('Features');  
ylabel('Features');  
title('Pearson Correlation Matrix of Features with Flood Probability');

% 后面根据需要补充 ... ...

方法2:使用斯皮尔曼等级相关系数(如果数据不服从正态分布)

因为数据可能不服从正态分布,所以需要计算每个特征与洪水概率之间的斯皮尔曼等级相关系数,将关系可视化便于识别每个系数之间特征。

% 加载数据  
data_train = readtable('train.xlsx');  
  
% 提取特征和标签(转换为数值数组因为这样便于计算斯皮尔曼相关系数)  
feature_matrix = table2array(data_train{:, 2:end-1});  
label_vector = data_train.洪水概率;  
  
% 计算斯皮尔曼等级相关系数  
rho_matrix = corr(feature_matrix, label_vector, 'Type', 'Spearman');  
  
% 注意:斯皮尔曼相关系数是特征与标签之间的相关系数,不是特征之间的。  
% 因此,又可能需要单独计算每个特征的相关系数,然后可视化它们。  
% 但为了与皮尔逊相关系数矩阵相似,所以可以考虑创建一个对角线为斯皮尔曼系数的矩阵  
num_features = size(feature_matrix, 2);  
corr_matrix_spearman = zeros(num_features);  
for i = 1:num_features  
    corr_matrix_spearman(i, i) = rho_matrix(i);  
end  
  
% 可视化斯皮尔曼等级相关系数(这里也可以选择使用条形图,都是可以的)
figure;  
bar(diag(corr_matrix_spearman));  
set(gca, 'XTick', 1:num_features);  
set(gca, 'XTickLabel', string(data_train.Properties.VariableNames(2:end-1)));  
xlabel('Features');  
ylabel('Spearman Rank Correlation Coefficient');  
title('Spearman Rank Correlation Coefficients of Features with Flood Probability');  
  
% 后面根据需要补充 ... ...

第二小问:聚类分析并建立预警评价模型

问题2原话“将附件train.csv 中洪水发生的概率聚类成不同类别,分析具有高、中、低风险的洪水事件的指标特征”

问题2已经启发我们了要聚类,所以这题突破口其实为找一个合适的聚类算法来解决该问题,我选择使用K-means聚类算法,因为洪水灾害数据包含大量的连续型特征,适合使用基于度量的聚类算法进行分析,同时在使用K-means聚类算法的同时我们可以结合方差检验来观测每个数据之间的差异。

聚类分析代码

% 加载数据  
data_train = readtable('train.xlsx');  
  
% 提取特征数据
features = data_train{:, 2:end-1};  % 假设第一列是'id',最后一列是'洪水概率'  
  
% 将table转换为数组
feature_matrix = table2array(features);  
  
% 使用肘部法则确定最佳的簇数量k  
wss = zeros(1, 10); % 只是假设1-10
for i = 1:10  
    [idx, C] = kmeans(feature_matrix, i);  
    wss(i) = sum(sum(bsxfun(@minus, feature_matrix, C(idx,:)).^2));  
end  
  
% 绘制肘部曲线图来确定k的最佳值  
figure;  
plot(1:10, wss);  
title('Elbow Method for Optimal Number of Clusters');  
xlabel('Number of Clusters');  
ylabel('Within-Cluster Sum of Squares');  
  
% 根据肘部法则,假设选择k=3作为最佳簇数量  
k = 3;  
  
% 使用选定的k值进行K-means聚类  
[idx, C] = kmeans(feature_matrix, k);  
  
% 将聚类结果添加到原始数据表中
data_train.ClusterIndex = idx;  
  
% 分析每个簇的特征 
% 计算每个簇的洪水概率平均值  
cluster_stats = groupsummary(data_train, 'ClusterIndex', {'mean(洪水概率)'});  
disp(cluster_stats);  
  
figure;  
for i = 1:k  
    subplot(k, 1, i);  
    histogram(feature_matrix(idx == i, 1)); % 假设我们查看第一个特征的分布  
    title(sprintf('Feature 1 Distribution in Cluster %d', i));  
end  

预警评价模型代码

% 直接使用前面代码中的结果,即data_train表包含了ClusterIndex列  
  
% 提取特征和簇索引  
features = data_train{:, 2:end-1};  % 假设特征从第二列开始,最后一列是洪水概率  
clusterIndex = data_train.ClusterIndex;  
  

% 直接定义一个权重向量作为示例  
featureWeights = [0.1, 0.05, 0.15, 0.1, 0.08, 0.12, 0.05, 0.07, 0.06, 0.09, ...  
                  0.04, 0.03, 0.06, 0.07, 0.02, 0.05, 0.04, 0.08, 0.03, 0.01];  
  
feature_matrix = table2array(features);  
warningScores = zeros(size(clusterIndex));  
for i = 1:max(clusterIndex)  
    currentCluster = feature_matrix(clusterIndex == i, :);   
    % 这个地方仅作为示例,实际中可能需要根据簇特征进行调整  
    for j = 1:size(currentCluster, 1)  
        warningScores(clusterIndex == i, j) = sum(currentCluster(j, :) .* featureWeights);  
    end  
end  
  
% 将预警评分添加到原始数据表里面
data_train.WarningScore = warningScores;  
   
% 制定相应的预警措施和应急预案。  
% 比如说,按照设定阈值来划分风险等级  
lowRiskThreshold = 0.3;  
highRiskThreshold = 0.7;  
data_train.RiskLevel = cellstr(arrayfun(@(x) ...  
    if x < lowRiskThreshold  
        'Low'  
    elseif x > highRiskThreshold  
        'High'  
    else  
        'Medium'  
    end, warningScores, 'UniformOutput', false));  

head(data_train, 10)

第三小问:建立洪水发生概率的预测模型

第一小问的指标分析结果就可以用在这个部分了,我们可以选择合适的特征建立洪水发生概率的预测模型,可以使用机器学习算法,比如说逻辑回归、随机森林等模型进行训练,因为train.csv中数据量很大便于我们训练模型

随机森林模型

假设已经从train.csv中加载了数据

% 加载数据  
data_train = readtable('train.xlsx');  
  
% 提取特征和目标变量  
% 假设特征从第二列到倒数第二列,最后一列是洪水概率  
X_train = table2array(data_train{:, 2:end-1});  % 特征矩阵  
y_train = data_train.洪水概率;  % 目标变量(洪水概率)  
  
% 分割数据为训练集和验证集(这里为了简化,我们直接使用全部数据作为训练集)  
model = fitlm(X_train, y_train);  
disp(model);  
  
% 使用验证集评估模型性能  
% 假设有验证集 X_val 和 y_val  
% y_pred = predict(model, X_val);  
% % 计算性能指标,如RMSE、MAE等  
  
% 由于小问3还提到了如果仅用5个关键指标的情况,这里我们假设已经通过某种方法(如特征选择)  
% 确定了5个关键指标,并重新建立模型  
% 假设关键指标的索引为 [1, 3, 5, 7, 9]
key_indices = [1, 3, 5, 7, 9];  
X_train_reduced = X_train(:, key_indices);  
  
% 使用关键指标重新建立线性回归模型  
model_reduced = fitlm(X_train_reduced, y_train);  
  
% 显示摘要  
disp(model_reduced);  
 
% 首先加载test数据  
data_test = readtable('test.xlsx');   
X_test = table2array(data_test{:, 2:end});  
y_pred_full = predict(model, X_test);  
y_pred_reduced = predict(model_reduced, X_test(:, key_indices));  
   
% 由于submit.xlsx文件中只有id列,需要创建一个新的表格来保存预测结果  
submit_table = table(data_test.id, y_pred_reduced, 'VariableNames', {'id', '洪水概率'});  
writetable(submit_table, 'submit.xlsx');

第四小问:预测test数据并分析结果分布

针对小问4的要求,需要基于前面建立的洪水发生概率的预测模型,对test.xlsx中的所有事件进行洪水概率的预测,并将预测结果填入submit.xlsx文件中,然后我们将绘制这些预测结果的直方图和折线图,并分析其分布是否服从正态分布。

% 加载test数据  
data_test = readtable('test.xlsx');  
  
% 这里直接使用那个模型进行预测  
  
% 提取特征
X_test = table2array(data_test{:, 2:end});  
  
% key_indices = [之前确定的关键指标索引];  
% X_test_reduced = X_test(:, key_indices);  
  
% 进行预测  
y_pred_full = predict(model, X_test);   
% y_pred_reduced = predict(model_reduced, X_test_reduced);  

% 创建包含id和预测洪水概率的表格  
submit_table = table(data_test.id, y_pred_full, 'VariableNames', {'id', '洪水概率'});  
  
% 写入Excel文件  
writetable(submit_table, 'submit.xlsx');

% 绘制直方图  
figure;  
histogram(y_pred_full, 'Normalization', 'pdf');  
title('预测洪水概率的直方图');  
xlabel('洪水概率');  
ylabel('频率密度');  
  
% 为了判断是否符合正态分布,可以叠加正态分布曲线
[mu, sigma] = normfit(y_pred_full);  
x = linspace(min(y_pred_full), max(y_pred_full), 1000);  
hold on;  
plot(x, normpdf(x, mu, sigma), 'r-');  
legend('数据直方图', '正态分布拟合');  
  
% 绘制折线图(通常对于大量数据点的连续概率值,折线图不是最佳选择,但这里为了完整性提供)  
figure;  
plot(data_test.id, y_pred_full, '-o');  
title('预测洪水概率的折线图');  
xlabel('事件ID');  
ylabel('洪水概率');  

% 后面根据需要补充 ... ...

B题matlab详细代码

B题python详细代码

作者:李 忘 忧

物联沃分享整理
物联沃-IOTWORD物联网 » 第十四届APMCM亚太地区大学生数学建模竞赛B题解题思路深度解析

发表回复