收藏！常用算法及MATLAB代码

Posted 2021-04-20 数学建模andMATLAB

 

主成分分析

主成分分析通常是选出比原始变量个数少，能解释大部分资料中的变异的几个新变量，即所谓主成分，并用以解释资料的综合性指标。主成分分析实际上是一种降维方法。

clc

clear

%累积贡献率阈值

p_threshold = 0.9;

%s输入原始数据

data = [149.3,4.2,108.1;161.2,4.1,114.8

           171.5,6.1,123.2;175.5,3.1,126.9

           180.8,8.1,132.1;110.7,2.2,137.7

           202.1,2.1,146.0;212.4,5.6,154.1

           226.1,5.0,162.3;231.9,5.1,164.3

           239.0,0.7,167.6];

[m,n] = size(data);

%数据的标准化处理

normal_data =(data -  repmat(mean(data),m,1)) ./repmat(std(data),m,1);

%计算协方差矩阵

sigama = cov(normal_data);

%计算协方差矩阵的特征值和特征向量

[V,lamadas] = eig(sigama);

%将特征值保存在一个行向量中

lamada = sum(lamadas);

%特征值进行降序排列，并记录特征值与特征向量的对应顺序

[order_lamada,index] = sort(-lamada);

order_lamada = -order_lamada;

%确定主成分个数

for i=1:length(lamada)

    P =sum(order_lamada(1:i)) / sum(order_lamada);%累计贡献率

    if P> p_threshold

       break

    end

end

num_pca = i;

%提取主成分的特征向量

V_main = V(:,index(1:i));

%计算主成分得分

new_score=normal_data * V_main;

%输出结果

disp('特征值、累计贡献率：')

order_lamada,P

disp('主成分个数与特征向量：')

num_pca

V_main

plot3(normal_data(:,1),normal_data(:,2),normal_data(:,3),'b*')

hold on

new_data = (V_main*V_main'*normal_data')';

plot3(new_data(:,1),new_data(:,2),new_data(:,3),'r--o')

xlabel('X');ylabel('Y');zlabel('Z');

title('原始数据和主成分得分对比')

legend('原始数据','主成分得分')

hold on

h =quiver3(0,0,0,V_main(1,1),V_main(2,1),V_main(3,1),2,'k','linewidth',1.5);

set(h,'maxheadsize',7);

hold on

h =quiver3(0,0,0,V_main(1,2),V_main(2,2),V_main(3,2),2,'k','linewidth',1.5);

set(h,'maxheadsize',4);

hold off

view(-27,10)

 

聚类分析

聚类分析又称群分析，是对多个样本（或指标）进行定量分类的一种多元统计分析方法。对样本进行分类称为Q型聚类分析，对指标进行分类称为R型聚类分析。

K-means聚类

clc

clear

close all

%输入原始数据

data =[0.697,0.460;0.774,0.376;0.634,0.264;0.608,0.318;0.556,0.215;0.403,0.237;

       0.481,0.149;0.437,0.211;0.666,0.091;0.243,0.267;0.245,0.057;0.343,0.099;

       0.639 0.161;0.657,0.198;0.360,0.370;0.593,0.042;0.719,0.103;0.359,0.188;

       0.339,0.241;0.282,0.257;0.748,0.232;0.714,0.346;0.483,0.312;0.478,0.437;

       0.525,0.369;0.751,0.489;0.532,0.472;0.473,0.376;0.725,0.445;0.446,0.459;];

%计算样本数量和特征数量

[num,num_feature] = size(data);

%聚类数量

k = 3;

%对每个特征数据按照极差法进行归一化

max_feature = max(data);

min_feature = min(data);

range_feature = max_feature - min_feature;

normal_data=(data - repmat(min_feature,num,1))./ repmat(range_feature,num,1);

%随机选择初始聚类中心

center = rand(k,num_feature);

%最大迭代次数

max_items = 10; count = 1;

while ( count <= max_items )

       %计算每个样本点到每个聚类中心的欧氏距离

       for  i = 1:k

           distance(i,:) = sum((normal_data - repmat(center(i,:),num,1))'.^2);

       end

       SSE = sum(min(distance))

       %判断每个样本点属于哪一类

       [~,belong_class] = min(distance);

       %更新聚类中心位置

       for i = 1:k

           center(i,:) = mean(normal_data(belong_class == i , :));

       end

       %记录每次迭代的聚类中心位置

       record_center(count,:) = reshape(center,1,k*num_feature);

       count = count + 1;

end

%聚类中心反归一化

center = repmat(min_feature,k,1) + center .* repmat(range_feature,k,1);

%判断聚类中心是否收敛

plot(record_center,'linewidth',2)

xlabel('迭代次数');

ylabel('聚类中心位置');

title('聚类中心位置变化趋势');

%描述不同样本的分类情况

figure

colormap = char(['r','b','k','g']);

fontmap = char(['*','^','o','p']);

for  i =1:k

   plot(data(belong_class == i,1),data(belong_class ==i,2),[colormap(i),fontmap(1)])

    holdon

   plot(center(i,1),center(i,2),[colormap(i),fontmap(2)])

end

xlabel('迭代次数'); ylabel('样本点'); title('样本点分类');

 

层次分析法

由专家比较两两指标的重要程度，主观确定权重的方法。

A =[1,1/5,1/2,1/3,1/2;5,1,3,4,3;2,1/3,1,2,1;3,1/4,1/2,1,2;2,1/3,1,1/2,1];

[n,~] = size(A);

%求出特征向量和特征值并求出最大特征值和它所对应的特征向量

[V,D] = eig(A);

[max_lamada , index ] = max(sum(D));

w = abs(V(:,index));

w = w ./ sum(w);

%以下是一致性检验

CI = (max_lamada - n ) / ( n - 1 );

RI = [0 0 0.52 0.89 1.12 1.26 1.36 1.41 1.461.49 1.52 1.54 1.56 1.58 1.59 1.60 1.61 1.615 1.62 1.63];

CR = CI / RI(n);

if CR<0.10

   disp('此矩阵的一致性可以接受!');

   disp('CI=');disp(CI);

   disp('CR=');disp(CR);

   disp('权重向量：');

   disp(w);

   disp('准则层最大特征根：');

   disp(max_lamada);

else

   disp('此矩阵的一致性验证失败，请重新进行评分!');

end

 

熵值法

根据指标的离散程度，客观确定权重的方法。

%实现用熵值法求各权重

clc

clear

%x为原始数据矩阵

data = [149.3,4.2,108.1;161.2,4.1,114.8

           171.5,6.1,123.2;175.5,3.1,126.9

           180.8,8.1,132.1;110.7,2.2,137.7

           202.1,2.1,146.0;212.4,5.6,154.1

           226.1,5.0,162.3;231.9,5.1,164.3

           239.0,0.7,167.6];

[m,n] = size(data); % n个样本, m个指标

%flag指示向量，1表示极大型指标，0表示极小型指标

flag = [1,1,0];

%%数据的归一化处理

max_feature = max(data);

min_feature = min(data);

range_feature = max_feature - min_feature;

normal_data=0.0001 + (data -repmat(min_feature,m,1)) ./ repmat(range_feature,m,1);

for i =1:n

    ifflag(i) == 0

        %极小型指标极大化处理

       normal_data(:,i) = 1- normal_data(:,i);

    end

end

%计算第j个指标下，第i个样本占该指标的比重p(i,j)

p = normal_data ./repmat(sum(normal_data),m,1);

%计算第j个指标的熵值e(j)

k = 1/log(n);

for j = 1:n

    e(j)= -k*sum(p(:,j).*log(p(:,j)));

end

d = ones(1,n) - e; %计算信息熵冗余度

w = real(d ./ sum(d)); %求权值w

遗传算法

遗传算法是一类借鉴生物界的进化规律(适者生存，优胜劣汰遗传机制)演化而来的随机化搜索方法。借鉴生物进化论，遗传算法将要解决的问题模拟成一个生物进化的过程，通过复制、交叉、突变等操作产生下一代的解，并逐步淘汰掉适应度函数值低的解，增加适应度函数值高的解。这样进化N代后就很有可能会进化出适应度函数值很高的个体。

 

clc

clear

close

population_size=100;%种群大小

chrom_length=10;%二进制编码长度

pc = 0.8;%交叉概率

pm = 0.01;%变异概率

%初始种群

population =round(rand(population_size,chrom_length));

max_iteration = 20;

for k = 1:max_iteration

    %计算适应度值（函数值）

   fitvalue = fitness(population);

    %选择操作

   new_population = selection(population,fitvalue);

    %交叉操作

   new_population = crossover(new_population,pc);

    %变异操作

new_population= mutation(new_population,pm);

     %更新种群

    population = new_population;

    %寻找最优解

    [max_fit,index] = max(fitvalue);

    current_best_pop = population(index,:);

    for i = 1:length(current_best_pop)

        current_best_pop_decimal(i) =2.^(length(current_best_pop)-i).*current_best_pop(i);

    end

    current_best_x(k) =sum(current_best_pop_decimal)*10/1023;

    if k > 1

        [current_best_y(k),idx] =max([current_best_y,max_fit]);

        current_best_x(k) =current_best_x(idx);

    else

        current_best_y(k) = max_fit;

    end

end

plot(current_best_y)

function F =fitness( population )

    [~,num2]=size(population);

    for i = 1:num2

        population_decimal(:,i) = 2.^(num2-i).*population(:,i);

    end

    %对行求和，得到列向量

    population_decimal =sum(population_decimal,2);

    x = population_decimal*10/1023;

    F = 10*sin(5*x)+7*abs(x-5)+10;

end

functionnew_population = selection(population,fitvalue)

    %构造轮盘

    [num1,~] = size(population);

    fitness_sum = sum(fitvalue);

    %计算累积概率

    p_fitvalue = cumsum(fitvalue /fitness_sum);

    ms = rand(num1,1);%从小到大排列

    for i= 1:num1

        selected = sum(p_fitvalue <=ms(i))+1;

       new_population(i,:)=population(selected,:);

    end

end

functionnew_population = crossover(population,pc)

    [num1,num2] = size(population);

    for i= 1:2:num1-1

        if(rand<pc)

            rand_int = round(rand*num2);

            new_population(i,:) =[population(i,1:rand_int),population(i+1,rand_int+1:num2)];

            new_population(i+1,:) =[population(i+1,1:rand_int),population(i,rand_int+1:num2)];

        else

            new_population(i,:) =population(i,:);

            new_population(i+1,:) =population(i+1,:);

        end

    end

end

functionnew_population = mutation(population,pm)

    [num1,num2] = size(population);

    new_population = population;

    for i = 1:num1

        for j = 1:num2

            if(rand < pm)

                if new_population(i,j) == 0

                    new_population(i,j) = 1;

                else

                    new_population(i,j) = 0;

                end

            end

        end

    end

end

 

粒子群算法

粒子群优化算法就是基于鸟群在觅食过程中的这种个体和整体的行为规则的启发，提出的一种模拟鸟群行为规则的智能优化算法。

clc

clear

close all

E=0.000001;

maxnum=800;%最大迭代次数

narvs=2;%目标函数的自变量个数

particlesize=50;%粒子群规模

c1=2;%每个粒子的个体学习因子，加速度常数

c2=2;%每个粒子的社会学习因子，加速度常数

w=0.6;%惯性因子

vmax=5;%粒子的最大飞翔速度

v=2*rand(particlesize,narvs);%粒子飞翔速度

x=-300+600*rand(particlesize,narvs);%粒子所在位置

%定义适应度函数

fitness=inline('(x(1)^2+x(2)^2)/10000','x');

for i=1:particlesize

       f(i)=fitness(x(i,:));      

end

personalbest_x=x;

personalbest_faval=f;

[globalbest_faval,i]=min(personalbest_faval);

globalbest_x=personalbest_x(i,:);

k=1;

while (k<=maxnum)

       fori=1:particlesize

              f(i)=fitness(x(i,:));

              iff(i)<personalbest_faval(i)

                     personalbest_faval(i)=f(i);

                     personalbest_x(i,:)=x(i,:);

              end

       end

       [globalbest_faval,i]=min(personalbest_faval);

       globalbest_x=personalbest_x(i,:);

       fori=1:particlesize

              v(i,:)=w*v(i,:)+c1*rand*(personalbest_x(i,:)-x(i,:))...

                     +c2*rand*(globalbest_x-x(i,:));

              forj=1:narvs

                     ifv(i,j)>vmax

                            v(i,j)=vmax;

                     elseifv(i,j)<-vmax

                            v(i,j)=-vmax;

                     end       

              end

              x(i,:)=x(i,:)+v(i,:);

    end

   ff(k)=globalbest_faval;

    ifglobalbest_faval<E

       break

    end

       k=k+1;

end

xbest=globalbest_x;

plot(1:length(ff),ff)

 

回归分析

回归分析(regression analysis)是确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法。

       自变量X：0.10,0.11,0.12,0.13,0.14,0.15,0.16,0.17,0.18

       自变量Y：42.0,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0

x=0.1:0.01:0.18;

y=[42,41.5,45.0,45.5,45.0,47.5,49.0,55.0,50.0];

%根据散点图观察变量之间关系

figure

plot(x,y,'r*')

xlabel('X'),ylabel('Y')

x=[ones(9,1),x'];

[b,bint,r,rint,stats]=regress(y',x);

figure

b,bint,stats,rcoplot(r,rint)

 

时间序列预测

移动平均法是根据时间序列的推移，依次计算包含一定项数的时序平均数，以反映长期趋势的方法。当时间序列的数值由于受周期变动和不规则变动的影响，起伏较大，不易显示出发展趋势时，可用移动平均法，消除这些因素的影响，分析、预测序列的长期趋势。

简单移动平均法

clc

clear

y=[533.8,574.6,606.9,649.8,705.1,772.0,816.4,892.7,963.9,1015.1,1102.7];

m=length(y);

n=[4,5]; %n 为移动平均的项数

for i=1:length(n)

    %由于 n 的取值不同， yhat 的长度不一致，所以使用元胞数组

    forj=1:m-n(i)+1

       yhat{i}(j)=sum(y(j:j+n(i)-1))/n(i);

    end

   y12(i)=yhat{i}(end);

   s(i)=sqrt(mean((y(n(i)+1:m)-yhat{i}(1:end-1)).^2));

end

y12,s

 

趋势移动平均法

简单移动平均法和加权移动平均法，在时间序列没有明显的趋势变动时，能够准确反映实际情况。但当时间序列出现直线增加或减少的变动趋势时，用简单移动平均法和加权移动平均法来预测就会出现滞后偏差。因此，需要通过二次移动平均法进行修正。

clc

clear

y =[676,825,774,716,940,1159,1384,1524,1668,1688,1958,

2031,2234,2566,2820,3006,3093,3277,3514,3770,4107];

m1=length(y);

n=6; %n 为移动平均的项数

for i=1:m1-n+1

   yhat1(i)=sum(y(i:i+n-1))/n;

end

m2=length(yhat1);

for i=1:m2-n+1

   yhat2(i)=sum(yhat1(i:i+n-1))/n;

end

plot(y,'*')

a21=2*yhat1(end)-yhat2(end);

b21=2*(yhat1(end)-yhat2(end))/(n-1);

y1986=a21+b21;

y1987=a21+2*b21;

一次指数平滑法

clc

clear

yt = [50,52,47,51,49,48,51,40,48,52,51,59]';

n=length(yt);

alpha=[0.2 0.5 0.8];

m=length(alpha);

yhat(1,1:m)=(yt(1)+yt(2))/2;

for i=2:n

    yhat(i,:)=alpha*yt(i-1)+(1-alpha).*yhat(i-1,:);

end

err=sqrt(mean((repmat(yt,1,m)-yhat).^2));

yhat1988=alpha*yt(n)+(1-alpha).*yhat(n,:);

二次指数平滑法

clc

clear

yt =[676,825,774,716,940,1159,1384,1524,1668,1688,

1958,2031,2234,2566,2820,3006,3093,3277,3514,3770,4107]';

n=length(yt);

alpha=0.3;

st1(1)=yt(1);

st2(1)=yt(1);

for i=2:n

   st1(i)=alpha*yt(i)+(1-alpha)*st1(i-1);

   st2(i)=alpha*st1(i)+(1-alpha)*st2(i-1);

end

a=2*st1-st2;

b=alpha/(1-alpha)*(st1-st2);

yhat=a+b;

 

灰色预测

通过鉴别系统因素之间发展趋势的相异程度，即进行关联分析，并对原始数据进行生成处理来寻找系统变动的规律，生成有较强规律性的数据系列，然后建立相应的灰色预测模型，对具有灰色系统特征的社会、经济等现象进行预测。

clc

clear

x0=[174 179 183 189 207 234 220.5 255 270285];

%输入原序列

n=length(x0);

%级比检验模型可行性

rank_ratio=zeros(1,n-1);

for i=2:n

   rank_ratio(i-1)=x0(i)/x0(i-1);

end

if (all(rank_ratio>=

exp(-2/(n+1)))&&all((rank_ratio<=exp(2/(n+1)))))

   fprintf('\n')

   disp('此序列适合GN（1，1）模型')

   fprintf('\n')

else

   disp('此序列不适合GN（1，1）模型')

end

x=cumsum(x0); %原序列做累加处理

B=zeros(1,n-1);

%产生均值序列

for i=1:n-1

   B(1,i)=(x(1,i)+x(1,i+1))/2;

end

B=[-B;ones(1,n-1)];

B=B';

Y=x0;

Y(1)=[];

Y=Y';

%最小二乘法求参数

a=inv(B'*B)*B'*Y;

F=zeros(1,20);F(1)=x0(1);

for i=2:n+10

   F(i)=(F(1)-a(2,1)/a(1,1))*exp(-a(1,1)*(i-1))+a(2,1)/a(1,1);

end                  %计算预测后累加序列

G=zeros(1,20);G(1,1)=x0(1,1);

for i=2:n+10

  G(i)=F(i)-F(i-1);

end                 %计算预测值

disp('输出预测值');

G

%相对误差检验

D=G(1,1:10);

res=abs(D-x0);%相对误差序列

rel_rate=res./x0;

rel_rate=mean(rel_rate);%计算相对误差

if rel_rate<=0.01

         disp('相对误差检验结果很好');

   elseif rel_rate<=0.05

        disp('相对误差检验结果合格');

   elseif rel_rate<=0.1

       disp('相对误差检验结果勉强合格');

end

fprintf('\n相对误差rel_rate=%4f\n\n',rel_rate);

%后验差检验

x0_mean=mean(x0);%原始序列均值

x0_std=std2(x0);%原始序列均方差

res_mean=mean(res);%残差序列均值

res_std=std2(res);%残差序列均方差

c=res_std/x0_std;%计算后验差比

s0=0.6745*x0_std;

e=abs(res-res_mean);

P=length(find(e<s0))/length(e);%计算小误差概率

if (c<=0.35&&P>=0.95)

   disp('后验差检验结果好');

end

if(c>0.35&&c<=0.5&&P>=0.8&&P<=0.95)

   disp('后验差检验合格');

end

if(c>0.5&&c<=0.65&&P>=0.7&&P<0.8)

   disp('后验差检验勉强合格');

end

fprintf('\n后验差比c=%4f\n',c);

fprintf('\n小误差概率P=%4f\n\n',P);

%关联度检验

abs_ss=D-x0;%绝对误差序列

k=(min(abs_ss)+0.5*max(abs_ss))./(abs_ss+0.5*max(abs_ss));

%计算关联系数

kk=mean(k);%计算平均 关联系数

if (abs(kk)>0.6)

   disp('关联度检验结果良好');

else

   disp('关联度检验不合格');

end

fprintf('\n关联度为%f\n\n',kk);

%残差检验

f=abs(D-x0)./x0;%残差序列

if (all(f<0.1))

        disp('残差检验结果好');

   elseif(all(f<0.2))

        disp('残差检验结果合格');

end

P0=1-sum(f)/(n-1);

fprintf('\n本次建模精度P0=%f\n',P0);