21-B题 学习笔记
这题今天没太做明白，这只是部分代码。
主要反思如下：应该观察全部21副图，或者看同一指标21条线的那两张关键图，然后发现由先上升后下降的趋势，因而应该选用二次函数拟合，即用一元二次函数建立回归模型利用最小二乘法求解（polyfit即可），而不是用一次函数。
其次，对于多变量相互影响的第二题，应该对数据做标准化处理，消除量纲和数量级带来的影响，然后建立考虑了所有量的方程，利用stepwise逐步回归进行求解更为合理,stepwise可以较好的解决单变量关于多变量做回归分析的问题。
新做法已经写完了，但不打算更新在这篇文章中，直接更新到《回归分析专题》中去了，下面这个完全不成熟的做法已经被取代废弃掉了
%导入
% 设置导入选项并导入数据
opts2 = spreadsheetImportOptions("NumVariables", 10);
% 指定工作表和范围
opts2.Sheet = "性能数据表";
opts2.DataRange = "A2:J115";
% 指定列名称和类型
opts2.VariableNames = ["x_______", "x_____", "x__", "x________", "x________1", "C4________", "x________2", "x___4_12_________", "x_________________", "x____________"];
opts2.VariableTypes = ["char", "char", "double", "double", "double", "double", "double", "double", "double", "double"];
% 指定变量属性
opts2 = setvaropts(opts2, ["x_______", "x_____"], "WhitespaceRule", "preserve");
opts2 = setvaropts(opts2, ["x_______", "x_____"], "EmptyFieldRule", "auto");
% 导入数据
Maintable = readtable("/Users/Kai/Desktop/B 2/附件1.xlsx", opts2, "UseExcel", false);
% 转换为输出类型
Maintable = table2cell(Maintable);
numIdx = cellfun(@(x) ~isnan(str2double(x)), Maintable);
Maintable(numIdx) = cellfun(@(x) {str2double(x)}, Maintable(numIdx));
% 清除临时变量
clear numIdx
% 清除临时变量
clear opts2
% 显示结果
%Maintable
index=["温度","乙醇转化率(%)","乙烯选择性（%）","C4烯烃选择性(%)","乙醛选择性(%)", ...
"碳数为4-12脂肪醇选择性(%)","甲基苯甲醛和甲基苯甲醇选择性(%)","其他生成物的选择性(%)"];
catalyst=strings(21,1);
table=cell(21,1);
p=1;
for i=1:21
catalyst(i,1)=Maintable{p,2};
lp=p;
p=p+1;
while(p<=114&&isempty(Maintable{p,2}))p=p+1;end
B=zeros(p-lp,8);
for j=lp:p-1
for k=3:10
B(j-lp+1,k-2)=Maintable{j,k};
end
end
table{i}=B;
end
%导入end
%程序
totaldata=1;
p1=ones(21,4);p2=p1;
needdraw=false;
for i=1:21
%提取数据
A=table{i};
T=A(:,1);
S=size(T);
tno=S(1); %本组数据所含温度数
%提取end
%绘制原始数据
if(needdraw)
for j=2:8
data=A(:,j);
hold on;
scatter(T,data,"LineWidth",2,"Marker",'o');
axis auto;
title("第"+num2str(i)+"组数据的图像");
xlabel('温度');
ylabel('选择性（%）')
xticks(T(1):25:T(tno));
yticks(0:10:100);
axis([T(1),T(tno),0,100]);
end
%绘图最后处理
legend(index(2:8)+"(原始数据)");%绘制图例
hold off;
%绘图最后处理end
end
%绘制原始数据end
C2H6O=A(:,2); %乙醇转化率
C4=A(:,4); %C4选择性
getC4=(C2H6O/100).(C4/100)100; %C4收率
p1(1:2,i)=polyfit(T,C2H6O,1);
p2(1:2,i)=polyfit(T,C4,1);
p3(1:2,i)=polyfit(T,getC4,1);
end
p1(1,:)
sum(p1(1,:))/21
p2(1,:)
sum(p2(1,:))/21
p3(1,:)
sum(p3(1,:))/21
%程序end
%合理认为温度越高、乙醇转化率和C4烯烃选择性越大,收率也越大
%取出所有乙醇浓度1.68ml/min的数据
cata=[200 200 1 1.68; ...
200 200 2 1.68; ...
200 200 1 0.9; ...
200 200 0.5 1.68; ...
200 200 2 0.3; ...
200 200 5 1.68; ...
50 50 1 0.3; ...
50 50 1 0.9; ...
50 50 1 2.1; ...
50 50 5 2.1; ...
50 90 1 1.68; ...
50 50 1 1.68; ...
67 33 1 1.68; ...
33 67 1 1.68; ...
50 50 1 1.68; ...
100 100 1 1.68; ...
10 10 1 1.68; ...
25 25 1 1.68; ...
50 50 1 2.1; ...
75 75 1 1.68; ...
100 100 1 0.9; ...
];
chooseindex=[1 2 4 6 11 12 13 14 15 16 17 18 20];
cata(:,1:2)=[cata(:,1)./cata(:,2),cata(:,3)];cata(:,3:4)=[];
choosecata=cata(chooseindex,1:2);
choosep1=p1(:,chooseindex)';
choosep2=p2(:,chooseindex)';
catax=choosecata(:,1);
catay=choosecata(:,2);
choosep1k=choosep1(:,1);
choosep1b=choosep1(:,2);
choosep2k=choosep2(:,1);
choosep2b=choosep2(:,2);
finaldata1=[catax catay choosep1k]
finaldata2=[catax catay choosep1b]
finaldata3=[catax catay choosep2k]
finaldata4=[catax catay choosep2b]
% trainingData=finaldata
%
% inputTable = array2table(trainingData, 'VariableNames', {'column_1', 'column_2', 'column_3', 'column_4', 'column_5'});
%
% predictorNames = {'column_1', 'column_2', 'column_3'};
% predictors = inputTable(:, predictorNames);
% response = inputTable.column_4;
% isCategoricalPredictor = [false, false, false];
%
% % 训练回归模型
% % 以下代码指定所有模型选项并训练模型。
% regressionNeuralNetwork = fitrnet(...
% predictors, ...
% response, ...
% 'LayerSizes', [10 10 10], ...
% 'Activations', 'relu', ...
% 'Lambda', 0, ...
% 'IterationLimit', 1000, ...
% 'Standardize', true);
%
% % 使用预测函数创建结果结构体
% predictorExtractionFcn = @(x) array2table(x, 'VariableNames', predictorNames);
% neuralNetworkPredictFcn = @(x) predict(regressionNeuralNetwork, x);
% trainedModel.predictFcn = @(x) neuralNetworkPredictFcn(predictorExtractionFcn(x));
%
% % 向结果结构体中添加字段
% trainedModel.RegressionNeuralNetwork = regressionNeuralNetwork;
% trainedModel.About = '此结构体是从回归学习器 R2023a 导出的训练模型。';
% trainedModel.HowToPredict = sprintf('要对新预测变量列矩阵 X 进行预测，请使用: \n yfit = c.predictFcn(X) \n将 ''c'' 替换为作为此结构体的变量的名称，例如 ''trainedModel''。\n \nX 必须包含正好 3 个列，因为此模型是使用 3 个预测变量进行训练的。\nX 必须仅包含与训练数据具有完全相同的顺序和格式的\n预测变量列。不要包含响应列或未导入 App 的任何列。\n \n有关详细信息，请参阅 How to predict using an exported model。');
%
% % 提取预测变量和响应
% % 以下代码将数据处理为合适的形状以训练模型。
% %
% % 将输入转换为表
% inputTable = array2table(trainingData, 'VariableNames', {'column_1', 'column_2', 'column_3', 'column_4', 'column_5'});
%
% predictorNames = {'column_1', 'column_2', 'column_3'};
% predictors = inputTable(:, predictorNames);
% response = inputTable.column_4;
% isCategoricalPredictor = [false, false, false];
%
% % 执行交叉验证
% partitionedModel = crossval(trainedModel.RegressionNeuralNetwork, 'KFold', 5)
%
% % 计算验证预测
% validationPredictions = kfoldPredict(partitionedModel)
%
% % 计算验证 RMSE
% validationRMSE = sqrt(kfoldLoss(partitionedModel, 'LossFun', 'mse'))
%
% trainedModel.predictFcn
% %neuralNetworkPredictFcn(predictorExtractionFcn(x))
% % 定义拟合函数
% fun = @(coefficients, x) coefficients(1) + coefficients(2)x(:, 1) + coefficients(3)x(:, 2) + coefficients(4)x(:, 3);
%
% % 初始参数估计
% beta1 = [0, 0, 0, 0];
% beta2 = beta1;
% w1=100;w2=10;%有待严谨化
%
% Y1=w1choosep1(:,1)+w2choosep1(:,2);
% Y2=w1choosep2(:,1)+w2*choosep2(:,2);
%
% % 拟合非线性回归模型
% model1 = fitnlm(X, Y1, fun, beta1);
% model2 = fitnlm(X, Y1, fun, beta2);
%
% % 获取参数估计值
% coefficients1 = model1.Coefficients;
% coefficients2 = model2.Coefficients;
%
% C1=table2array(coefficients1(:,1))
% C2=table2array(coefficients2(:,1))
clc;clear;close all;
%% 设置导入选项并导入数据
opts = spreadsheetImportOptions("NumVariables", 8);
% 指定工作表和范围
opts.Sheet = "稳定性测试";
opts.DataRange = "A4:H10";
% 指定列名称和类型
opts.VariableNames = ["VarName1", "VarName2", "VarName3", "C4", "VarName5", "VarName6", "VarName7", "VarName8"];
opts.VariableTypes = ["double", "double", "double", "double", "double", "double", "double", "double"];
% 导入数据
table2 = readtable("/Users/Kai/Desktop/B 2/附件2.xlsx", opts, "UseExcel", false);
%% 转换为输出类型
table2 = table2array(table2);
%% 清除临时变量
clear opts
given_T=350;
t=table2(:,1);
index=["温度","乙醇转化率(%)","乙烯选择性（%）","C4烯烃选择性(%)","乙醛选择性(%)", ...
"碳数为4-12脂肪醇选择性(%)","甲基苯甲醛和甲基苯甲醇选择性(%)","其他生成物的选择性(%)"];
S=size(t);
tno=S(1);
hold on;
for j=2:8
data=table2(:,j);
scatter(t,data,"LineWidth",2,"Marker",'o');
end
pol=ones(5,8);
for j=2:8
data=table2(:,j);
pol(1:4,j)=polyfit(t,data,3);
x=20:0.05:280;
y=polyval(pol(1:4,j),x);
plot(x,y,"LineWidth",3);
end
xlabel('时间');
ylabel('百分比');
xticks(linspace(20,280,14));
yticks(0:10:100);
axis([20,300,0,100]);
legend([index(2:8)+"(原始数据)",index(2:8)+"(拟合曲线)"]);%绘制图例
hold off;
pol(1:4,2:8)
