2021年中国研究生数学建模竞赛D题参考思路
问题一难度较低，主要是找出对生物活性最具有显著影响的分子描述符，具有影响意味着该分子描述符与生物活性的相关性越高，因此可以构建相关性分析模型来分别计算各分子描述符与生物活性的关系，找出相关性最大的20个变量即可。目前有同学说用皮尔逊系数做的结果并不是很好，这是因为皮尔逊系数主要适用于单调线性变化的数据，而题目中的数据并没有明显的单调性质，因此得到的结果可能并不好，这时候建议大家可以试一下灰色关联分析模型，他是根据因素之间发展趋势的相似或相异程度，亦即“灰色关联度”，作为衡量因素间关联程度的一种方法。
灰色关联分析代码如下
％数据替换为自身研究数据，代码会自动去除量纲
clear;clc;
load data.mat;
r = size(data,1);
c = size(data,2);
%第一步，对变量进行预处理，消除量纲的影响
avg = repmat(mean(data),r,1);
data = data./avg;
%定义母序列和子序列
Y = data(:,1); %母序列
X = data(:,2:c); %子序列
Y2 = repmat(Y,1,c-1); %把母序列向右复制到c-1列
absXi_Y = abs(X-Y2)
a = min(min(absXi_Y)) %全局最小值
b = max(max(absXi_Y)) %全局最大值
ro = 0.5; %分辨系数取0.5
gamma = (a+ro*b)./(absXi_Y+ro*b) %计算子序列中各个指标与母序列的关联系数
disp("子序列中各个指标的灰色关联度分别为：");
ans = mean(gamma)针对问题二：该问题主要是构建生物活性与20个分子变量之间的逻辑关系，主要有两种思路，第一种是构建多元线性模型，及构建y=a1x1+a2x2+...+a20x20关系，然后求解a1,a2,...a20系数即可；第二种思路是构建神经网络模型，建立变量与生物活性之间的非线性拟合关系，这里最好别用BP神经网络，因为输入数据的参数较多，很容易陷入局部最优解，可采用遗传算法优化的BP神经网络模型或粒子群算法优化的神经网络模型等，建立上述逻辑关系后将test表中数据带入即可。
有同学说这个关系如何构建，能不能用回归分析的思路，回归分析肯定是用线性回归，非线性对于这么多的指标难度太大而且也不好把握，但如果线性回归其关系式又太过简单（毕竟医学领域用线性关系构建的可能性很小）因此回归算法在这里并不推荐。
还是主推神经网络之类的非线性拟合模型，对于未知关系的拟合不借助任何参考关系式，而是采用权值链接在一起，精度会更高一些，但注意过拟合和误差分析问题，并不局限于用BP神经网络，径向基网络，LM－BP神经网络等也可以试一下
%程序一：GA训练BP权值的主函数
function net=GABPNET(XX,YY)
%--------------------------------------------------------------------------
% GABPNET.m
% 使用遗传算法对BP网络权值阈值进行优化，再用BP算法训练网络
%--------------------------------------------------------------------------
%数据归一化预处理
nntwarn off
XX=[1:19;2:20;3:21;4:22]';
YY=[1:4];
XX=premnmx(XX);
YY=premnmx(YY);
YY
%创建网络
net=newff(minmax(XX),[19,25,1],{'tansig','tansig','purelin'},'trainlm');
%下面使用遗传算法对网络进行优化
P=XX;
T=YY;
R=size(P,1);
S2=size(T,1);
S1=25;%隐含层节点数
S=RS1+S1S2+S1+S2;%遗传算法编码长度
aa=ones(S,1)*[-1,1];
popu=50;%种群规模
save data2 XX YY % 是将 xx,yy 二个变数的数值存入 data2 这个MAT-file，
initPpp=initializega(popu,aa,'gabpEval');%初始化种群
gen=100;%遗传代数
%下面调用gaot工具箱，其中目标函数定义为gabpEval
[x,endPop,bPop,trace]=ga(aa,'gabpEval',[],initPpp,[1e-6 1 1],'maxGenTerm',gen,...
'normGeomSelect',[0.09],['arithXover'],[2],'nonUnifMutation',[2 gen 3]);
%绘收敛曲线图
figure(1)
plot(trace(:,1),1./trace(:,3),'r-');
hold on
plot(trace(:,1),1./trace(:,2),'b-');
xlabel('Generation');
ylabel('Sum-Squared Error');
figure(2)
plot(trace(:,1),trace(:,3),'r-');
hold on
plot(trace(:,1),trace(:,2),'b-');
xlabel('Generation');
ylabel('Fittness');
%下面将初步得到的权值矩阵赋给尚未开始训练的BP网络
[W1,B1,W2,B2,P,T,A1,A2,SE,val]=gadecod(x);
net.LW{2,1}=W1;
net.LW{3,2}=W2;
net.b{2,1}=B1;
net.b{3,1}=B2;
XX=P;
YY=T;
%设置训练参数
net.trainParam.show=1;
net.trainParam.lr=1;
net.trainParam.epochs=50;
net.trainParam.goal=0.001;
%训练网络
net=train(net,XX,YY);
%程序二：适应值函数
function [sol, val] = gabpEval(sol,options)
% val - the fittness of this individual
% sol - the individual, returned to allow for Lamarckian evolution
% options - [current_generation]
load data2
nntwarn off
XX=premnmx(XX);
YY=premnmx(YY);
P=XX;
T=YY;
R=size(P,1);
S2=size(T,1);
S1=25;%隐含层节点数
S=RS1+S1S2+S1+S2;%遗传算法编码长度
for i=1:S,
x(i)=sol(i);
end;
[W1, B1, W2, B2, P, T, A1, A2, SE, val]=gadecod(x);
%程序三：编解码函数
function [W1, B1, W2, B2, P, T, A1, A2, SE, val]=gadecod(x)
load data2
nntwarn off
XX=premnmx(XX);
YY=premnmx(YY);
P=XX;
T=YY;
R=size(P,1);
S2=size(T,1);
S1=25;%隐含层节点数
S=RS1+S1S2+S1+S2;%遗传算法编码长度
% 前RS1个编码为W1
for i=1:S1,
for k=1:R,
W1(i,k)=x(R(i-1)+k);
end
end
% 接着的S1S2个编码（即第RS1个后的编码）为W2
for i=1:S2,
for k=1:S1,
W2(i,k)=x(S1*(i-1)+k+RS1);
end
end
% 接着的S1个编码（即第RS1+S1S2个后的编码）为B1
for i=1:S1,
B1(i,1)=x((RS1+S1S2)+i);
end
% 接着的S2个编码（即第RS1+S1S2+S1个后的编码）为B2
for i=1:S2,
B2(i,1)=x((RS1+S1S2+S1)+i);
end
% 计算S1与S2层的输出
A1=tansig(W1P,B1);
A2=purelin(W2*A1,B2);
% 计算误差平方和
SE=sumsqr(T-A2);
val=1/SE; % 遗传算法的适应值
针对第三问：该部分与问题二思路大致相同，但由于ADMET数据都是由0和1组成，因此若采用神经网络模型时需要考虑隶属度问题以及负值问题，同时根据建模常识尽量不要用同一种算法求解全部关键问题，但可以是对算法的升级，因此这里也可以采用深度学习相关的分类算法来做，主流的深度学习方法包括：AlexNet、VGG-Net、ResNet等；同时在这里推荐使用其它的分类算法，例如SVM支持向量机、随机森林、adaboost、xgboost等。
对于第三问而言，主要是构建分类模型，这里有同学想的是构建五个不同的分类模型，比如综合用神经网络、SVM、、、等五个模型，其实完全没必要，因为该分类问题目标类别只有两种，但主流的分类方法都适合于多类别的划分，因此在设计模型时（尤其是神经网络模型，要注意隶属度的问题），但在这里重点推荐逻辑回归分类模型，目前Logistic回归主要就是用在二分类领域。它是一个非线性的回归模型， 其最大的好处恰恰是可以解决二元类问题，目前在很多行业，基本都是使用Logistic回归来预判目标是好还是坏，因为它还弥补了其他黑盒模型(SVM、 神经网络、随机森林等)不具解释性的缺点。
Logistic回归分类代码如下
分类案例：本文使用的例子是安德森鸢尾花卉数据集（iris），这个数据集中包含150个样本，对应着它的150行，每一行包含这个样本的4个特征（花萼长度、花萼宽度、花瓣长度、花瓣宽度）和样本的类别标签（0或1或2，它们分别代表不同的品种），所以数据集中存储的是一个尺寸为150×5的二维矩阵（数据集中部分样本的信息如下图所示）。由于我们此次只研究二分类问题，所以我们只选用数据集的前100个样本，即类别标签为“0”或“1”的样本。我们的任务是：建立一个逻辑回归分类器，这个分类器可以通过样本的4个特征值来预测样本的标签是“0”还是“1”。
iris数据集我已经上传到百度云，有需要的小伙伴可以自行下载。
链接：https://pan.baidu.com/s/15p8-JOoui8taai0tSADNyA
提取码：k0cs
以上数据集仅供参考，下面替换成自己的就行
% 读取数据集，变量data存储的是一个尺寸为150×5的矩阵
data = load('iris.data');
% 只取前100个样本的信息，即前100行，变量useful_data存储的是一个尺寸为100×5的矩阵
useful_data = data(1:100, :);
% 将特征与类别标签分开存放
% 特征存放在变量X中，X存储的是一个尺寸为100×4的矩阵
% 类别标签存放在变量y中，y存储的是一个尺寸为100×1的矩阵
X = useful_data(:, 1:4);
y = useful_data(:, 5);
% 变量m存储的是变量X的行数，在这里为100
% 变量n存储的是变量X的列数，在这里为4
[m,n] = size(X);
% 在变量X后加一列“1”，便于后面使用矩阵运算简化步骤
X = [X ones(m, 1)];
% 初始化模型参数β=(ω,b)为0
beta = zeros(n+1, 1);
% 设置梯度下降迭代次数为1500次
iteration = 1500;
% 设置学习率为0.01
alpha = 0.01;
% 开始循环，用梯度下降更新参数
for iter = 1 : iteration
z = X * beta;
h = 1 ./ (1 + exp(-z));
error = h - y;
graident = X' * error;
beta = beta - alpha / m * graident;
end
针对问题四：
问题4主要是探究哪些分子描述符对生物活性和ADMET效果更好，在这里我们可以按照两步走的方法，首先是利用问题一选择的20个变量和问题二构建的模型（在此以神经网络模型为例），构建完神经网络模型后，此时编写相关的程序如下：
fun = @(x) -predict(net, x);
x0 = [0,0,0,0]; %你自己设置初值，（有可能需要转置即x0 = [0,0,0,0]'）
[x, ymax] = fminsearch(fun, x0); x就是你要的位置
ymax = -ymax;%ymax就是极大值
此时就能找到能够使生物活性达到最大值时上述分子描述符的取值，然后将上述初值带入到问题三构建的分类模型中，看看ADMET分别取值是多少，若有三个性质较好则表示已经完成，若效果并不好，则返回第一步按照问题一中相关性高低原则逐渐剔除因子，指导既满足分子活性最大又使得ADMET中至少三个性质较好。
