【matlab 投影寻踪】基于PSO算法的最优投影方向优化

一投影寻踪算法

投影寻踪是处理和分析高维数据的一类统计方法，其基本思想是将高维数据投影到低维（1～3维）子空间上，寻找出反映原高维数据的结构或特征的投影，以达到研究和分析高维数据的目的。1974年，美国Stanford大学的Friedman和Tukey首次将该方法命名为Projection Pursuit，即投影寻踪。
投影寻踪（projection pursuit，简称PP）是国际统计界于70年代中期发展起来的一种新的、有价值的新技术，是统计学、应用数学和计算机技术的交叉学科。它是用来分析和处理高维观测数据，尤其是非正态非线性高维数据的一种新兴统计方法。它通过把高维数据投影到低维子空间上，寻找出能反映原高维数据的结构或特征的投影，达到研究分析高维数据的目的。它具有稳健性、抗干扰性和准确度高等优点，因而在许多领域得到广泛应用。

三粒子群算法

粒子群算法（Particle Swarm Optimization，PSO）是一种基于群体智能的优化算法，用于解决优化问题。它模拟群体中的个体（粒子）在解空间中的搜索过程，通过个体之间的信息交流和合作来寻找最优解。

粒子群算法的基本思想是，将解空间看作是一个多维空间，个体（粒子）在解空间中移动，并通过个体之间的信息交流和自我学习来改进自己的搜索策略。每个个体都有自己的位置和速度，并根据自己的经验（个体最优解）和群体的经验（全局最优解）来调整自己的移动方向和速度。

在粒子群算法中，每个粒子都有一个位置向量和一个速度向量。粒子根据自己的速度向量更新自己的位置，并根据自己的位置计算自己的适应度值。每个粒子都保留了自己的最优位置和适应度值，同时也保存了群体中的全局最优位置和适应度值。粒子根据自己的最优位置和全局最优位置来更新自己的速度，从而改变自己的移动方向和速度。

具体来说，粒子的速度更新公式为：
v_i(t+1) = w * v_i(t) + c1 * r1 * (pbest_i - x_i(t)) + c2 * r2 * (gbest - x_i(t))

其中，v_i(t)为粒子i在时刻t的速度，w为惯性权重，c1和c2为加速因子，r1和r2为随机数。pbest_i为粒子i的最优位置，x_i(t)为粒子i在时刻t的位置，gbest为全局最优位置。

粒子的位置更新公式为：
x_i(t+1) = x_i(t) + v_i(t+1)

粒子群算法的搜索过程一般包括以下步骤：

初始化粒子的位置和速度。
计算每个粒子的适应度值，并更新每个粒子的最优位置和全局最优位置。
根据粒子的最优位置和全局最优位置，更新每个粒子的速度和位置。
重复步骤2和步骤3，直到达到停止条件。

粒子群算法的优点是简单易实现，并且能够在较短的时间内找到较优解。然而，它也存在一些缺点，比如容易陷入局部最优解、对参数敏感等。因此，在应用粒子群算法时需要根据实际情况进行参数的调整和优化。

三代码分享

clc
close all
clear all
dbstop if error
%%
filename='一二三四五线城市量化分析数据.xlsx';
% sheetname='2000-2013';
sheetname='2014-2016';[data,txt]=xlsread(filename,sheetname);
%% 读取数据
% xij=data(:,9:end-2);
xij=data(:,9:end);namelist=txt(3:end,2);
x_ij_norm=(xij-min(min(xij)))./(max(max(xij))-min(min(xij)));
%% 求解最优投影系数a_best
[n,m]=size(x_ij_norm);
objfunction=@(a)objFunction(a,x_ij_norm);
%%求解
lb=zeros(m,1);
ub=ones(m,1);
pop=50;
Max_iter=100;
[a_best,Best_score,curve]=PSO(pop,Max_iter,lb,ub,m,objfunction);figureplot(abs(curve))
xlabel('迭代次数')
ylabel('目标函数指')
title('粒子群算法进化曲线')%% 计算落户门槛指数
Z=x_ij_norm*a_best';%% 显示结果figure
subplot(2,2,1)
bar(Z(1:30))
xticks(1:30)
xticklabels(namelist(1:30))
ylabel('落户门槛指数')
subplot(2,2,2)
bar(Z(1+30:60))
xticks(1:30)
xticklabels(namelist(1+30:60))
ylabel('落户门槛指数')
subplot(2,2,3)
bar(Z(1+60:90))
xticks(1:30)
xticklabels(namelist(1+60:90))
ylabel('落户门槛指数')
subplot(2,2,4)
bar(Z(1+90:120))
xticks(1:30)
xticklabels(namelist(1+90:120))
ylabel('落户门槛指数')
sgtitle(sheetname)

%% 粒子群算法
function [Best_pos,Best_score,curve]=PSO(pop,Max_iter,lb,ub,dim,fobj)
%% 参数设置
w = 0.9;      % 惯性因子
c1 = 2;       % 加速常数
c2 = 2;       % 加速常数
%速度范围设定
Vmax = 2;
Vmin = -2;
Dim = dim ;            % 维数
sizepop = pop;       % 粒子群规模
maxiter  = Max_iter;      % 最大迭代次数
%%
if(max(size(ub)) == 1)ub = ub.*ones(1,dim);lb = lb.*ones(1,dim);
end
fun=fobj;
%% 粒子群初始化
% Range = ones(sizepop,1)*(ub-lb);
% pop = rand(sizepop,Dim).*Range + ones(sizepop,1)*lb;    % 初始化粒子群
pop=init(sizepop,lb,ub);
V = rand(sizepop,Dim)*(Vmax-Vmin) + Vmin;                 % 初始化速度
fitness = zeros(sizepop,1);
for i=1:sizepopfitness(i,:) = fun(pop(i,:));                         % 粒子群的适应值
end%% 个体极值和群体极值
[bestf, bestindex]=min(fitness);
zbest=pop(bestindex,:);   % 全局最佳
gbest=pop;                % 个体最佳
fitnessgbest=fitness;              % 个体最佳适应值
fitnesszbest=bestf;               % 全局最佳适应值%% 迭代寻优
iter = 0;
while( (iter < maxiter ))for j=1:sizepop% 速度更新V(j,:) = w*V(j,:) + c1*rand*(gbest(j,:) - pop(j,:)) + c2*rand*(zbest - pop(j,:));if V(j,:)>VmaxV(j,:)=Vmax;endif V(j,:)<VminV(j,:)=Vmin;end% 位置更新pop(j,:)=pop(j,:)+V(j,:);for k=1:Dimif pop(j,k)>ub(k)pop(j,k)=ub(k);endif pop(j,k)<lb(k)pop(j,k)=lb(k);endendpop(j,:)=pop(j,:)./sum(pop(j,:));% 适应值fitness(j,:) =fun(pop(j,:));% 个体最优更新if fitness(j) < fitnessgbest(j)gbest(j,:) = pop(j,:);fitnessgbest(j) = fitness(j);end% 群体最优更新if fitness(j) < fitnesszbestzbest = pop(j,:);fitnesszbest = fitness(j);endenditer = iter+1;                      % 迭代次数更新fprintf('当前迭代次数%i\n',iter)curve(iter) = fitnesszbest;
end
%% 绘图
Best_pos = zbest;
Best_score = fitnesszbest;