一 投影寻踪算法
投影寻踪是处理和分析高维数据的一类统计方法，其基本思想是将高维数据投影到低维（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);
figure
plot(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:sizepop
fitness(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,:)>Vmax
V(j,:)=Vmax;
end
if V(j,:)<Vmin
V(j,:)=Vmin;
end
% 位置更新
pop(j,:)=pop(j,:)+V(j,:);
for k=1:Dim
if pop(j,k)>ub(k)
pop(j,k)=ub(k);
end
if pop(j,k)<lb(k)
pop(j,k)=lb(k);
end
end
pop(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) < fitnesszbest
zbest = pop(j,:);
fitnesszbest = fitness(j);
end
end
iter = iter+1; % 迭代次数更新
fprintf('当前迭代次数%i\n',iter)
curve(iter) = fitnesszbest;
end
%% 绘图
Best_pos = zbest;
Best_score = fitnesszbest;
四 结果演示
