美赛整理之投影寻踪模型及其求解

1.模型的简介和应用

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

2.基本步骤：

​ 设第$i$个样本的第$j$个指标为$x_{ij}(i=1,2..n,j=1,2..m)$，$n$为样本个数，$m$为指标个数。投影寻踪模型建立的投影寻踪聚类模型的步骤如下。

​ 1.样本数据归一化：
$$x_{ij}^{‘}=\frac{x_{ij}-x_{j,min}}{x_{j,max}-x_{j,min}}$$
​ 2.线性投影：

​ 从不同的角度去观察数据，寻找最能挖掘数据大的最优投影方向。随机抽取若干初始投影方向$a(a_1,a_2,...,a_m)$,计算其投影指标的大小。确定最大指标的投影的解为其最佳投影方向：

​ 样本$i$在一维空间的投影方向为$(a_1,a_2,...a_m)$上的投影特征值定义为：
$$z_i=\sum _{j=1}^m{a}_j{x}_{ij}$$
​ 3.寻找目标函数：

​ 要求$z_i$的为:局部投影点尽可能集中，形成许多的点团。整体上尽可能散开。

​ 定义目标函数为$Q(a)$,类间距离为$s(a)$,类内密度为$d(a)$。那么：
$$Q(a)=s(a)d(a)$$
​ 类间距离用样本序列的投影特征值方差计算：
$$\begin{gathered} s(a)=\sqrt{\frac{1}{n}\sum _{i=1}^n(z_i-\frac{1}{n}\sum _{i=1}^n{z}_i{)}^2} \\ i=1,2...n \end{gathered}$$
​ 投影特征值间的距离为$r_{ik}=|z_i-z_k|(i,k=1,2,...n)$，那么：
$$d(a)=\sum _{i=1}^n\sum _{k=1}^n(R-r_{ik})f(R-r_{ik})$$
​ 其中的$f(t)$为以阶跃信号:
$$f(t)=\{\begin{array}{l}0,t<0\\ 1,t\ge 0\end{array}$$
​ $R$为估计局部散点密度的窗宽参数，按宽度内至少包括一个散点的原则，其值与样本数据结构有关：
$$\begin{gathered} max(r_{ik})<R<2m \\ i,k=1,2,...n \end{gathered}$$
​ 4.优化投影方向：

​ 将寻找最优投影方向问题转化为下列优化问题：
$$\begin{gathered} maxQ(a)=s(a)d(a) \\ \sum _{j=1}^m{a}_j^2=1 \end{gathered}$$
​ 可用遗传算法优化。

​ 5.综合评价聚类分析：

​ 以$z_i$的差异水平对样本群进行聚类分析。根据最优投影方向，便可计算反应各评价指标综合信息的投影特征值$z_i$的差异水平。求出最优投影系数$a=a(a_1,a_2,...a_m)$。

3.遗传算法求解模型的优化问题

1.用$Matlab$的$ga$工具箱求解：

​ 调用$Matlab$的遗传算法工具箱。
$$[x,fval]=ga(@fitness,nvars,A,b,Aeq,beq,LB,UB,@nonlcon,options);$$
​ 其中$x$为决策向量，$fitness$为适应度函数$fitness(x)$，$nvars$：代表决策向量$x$的个数，其中$A,b,Aeq,beq,LB,UB$满足以下关系：
$$\begin{gathered} Ax\le b \\ Aeq\times x\le beq \\ LB\le x\le UB \end{gathered}$$
​ 其中$[f,g]=nonlcon(x)$，且$f$为非线性不等式约束，$g$为非线性等式约束:
$$\begin{gathered} f(x)\le 0 \\ g(x)=0 \end{gathered}$$
​ 而很多时候其实我们如果有非线性等式约束的话是很难直接找到合适的解的，或者说要花很长时间。这是我们可以设定一个阈值$\delta$,使得：
$$0\le g(x)\le \delta$$
​ 来大概满足我们的求解要求。那么此时就可以把等式约束化为不等式约束：
$$f=\left(\begin{array}{c}-g(x)\\ g(x)-\delta \end{array}\right)$$
​ 通过这种手段尽可能搜索最优解。相关的matlab代码如下：
归一化代码：

function [z,xd] = normalizef( x,a )
%UNTITLED6 此处显示有关此函数的摘要
%   此处显示详细说明
% [n1,m1]= size(a);%默认是一行m1列
 [n,m] = size(x);
x_min = zeros(1,m);
x_max = zeros(1,m);
z = zeros(1,n);
 for j = 1:m
    x_min(j) = min(x(:,j));
    x_max(j) = max(x(:,j));
end
for i = 1:n
    for j = 1:m
        xd(i,j) = (x(i,j)-x_min(j))/(x_max(j)-x_min(j));
    end
end
for i = 1:n
    sum = 0;
    for j = 1:m
        sum = sum + a(j)*xd(i,j);
    end
    z(i) = sum;
end

end

求$S(a)$的代码：

function S_a = S(z)
%UNTITLED2 此处显示有关此函数的摘要
%   此求S(a)
[m,n] = size(z);%默认为是一行n列
z_mean = sum(z)/n;
S_a = sqrt( 1/n*(sum((z-z_mean).^2)) );

end

求$D(a)$的代码：

function da = D(z)
%UNTITLED4 此处显示有关此函数的摘要
%   求d(a),我们假定R = r_m+3;
[m,n] = size(z);%默认为是一行n列
sum = 0;
for i = 1:n
    for k = 1:n
        r(i,k) = abs(z(i)-z(k));
        % sum = sum+(R-r(i,k))*fun(R-r(i,k));
    end
end

R = max(max(r));
R = R+3;

for i = 1:n
    for k = 1:n
       sum = sum+(R-r(i,k))*fun(R-r(i,k));
    end
end

da = sum;
end

函数$f(x)$的代码：

function f = fun( t )
%UNTITLED5 此处显示有关此函数的摘要
%   此处显示详细说明
if t>=0
    f =1;
else
    f = 0;
end
end

适应度函数代码：

function Q = fitness(a)
%UNTITLED8 此处显示有关此函数的摘要
%   用遗传算法求最优的 a (适应度函数)
x = xlsread('数据.xlsx','sheet1','B2：N9');
[z,xd] = normalizef( x,a );
 Q = S(z)*D(z);
 Q = -Q;
end

约束函数代码:

function [ f,g ] = yueshu( a )
%UNTITLED10 此处显示有关此函数的摘要
%   此处显示详细说明
%[n,m] = size(a);
delta = 0.1;
f = [1-sum(a.^2);sum(a.^2)-1-delta;];
g= [];
end

主函数代码:

%主函数用遗传算法求解最合适的投影角度
nvars = 13;
a = sqrt(1/13)*ones(1,nvars);
options = gaoptimset('PopulationSize',50,'Generations',20,'PlotFcns',@gaplotbestf);
[a fval] = ga(@fitness,nvars,[],[],[],[],[],[],@yueshu,options);
fval = -fval;
disp('最终的a是：');
disp(a);
disp('目标函数值是：');
disp(fval);

结果是:

Optimization terminated: average change in the fitness value less than options.FunctionTolerance
 and constraint violation is less than options.ConstraintTolerance.
最终的a是：
   -0.1534    0.2417   -0.0865    0.7742   -0.0963   -0.0297    0.2373    0.3051   -0.3150    0.3357   -0.0384    0.0405   -0.0820

经过验证得到:

sum(a.^2)

ans =

   1.0702

数据暂不给出

遗传算法工具箱的优化的确有点复杂，后面将直接用遗传算法来优化。

2.遗传算法求解：

如果用遗传算法求解的话最重要的就是要有一个函数能产生满足约束条件${\sum }_{i=1}^m{a}_i^2=1$的解$(a_1,a_2,...a_m)$。现在就将这样的满足约束条件的赋值函数给出：

function al = inital(nvars,N,step)
%nvars = 13;
%N = 10000;决定最后样本个数的次数
%step = 0.1;划分步长
al = [];
k = 0;
for alpha = 0:step:1
    for k=1:N
        a = zeros(1,nvars);
        k = k+1;
        for i = 1:(nvars-1)
            a(i) = rand*(1-alpha);
        end
        sum_n_1 = sum(a(1:(nvars-1)).^2);
        if sum_n_1>1
            continue;
        else
            a(nvars) = sqrt(1- sum_n_1);
            al = [al;a];
        end
    end
end
% [m,n] = size(al);
%  for j = 1:100
%      Q(j) = fitness(al(j,:));
%  end
%  plot(1:100,Q);
end

这个算法的主要思想是用不同比例的$[0,1]$的均匀随机数用来填充$(a_1,a_2,...a_{m-1})$,然后根据这$m-1$组数据求出$a_m$。这样就可以求出满足约束条件${\sum }_{i=1}^m{a}_i^2=1$的很多样本。然后由大数定律可以求出，最后实际满足要求的样本个数约为$\frac{N}{2step}$。

剩下的遗传算法的其他实现过程就不再赘述。
取前$100$组数据观察这些数据的适应度为:

3.$LINGO$求解

若将归一化后的数据存入$数据2.xlsx$的$data$区域后，调用$LINGO$求解：

model:
sets:
row/1..8/:z;
col/1..13/:a;
link(row,col):x;
same(row,row):r;
endsets
data:
x = @ole('D:\LINGO11\数据2.xlsx',data2);
enddata
max = s*d;
n = @size(row);
m = @size(col);
@for(row(i):z(i)=@sum(col(j):a(j)*x(i,j)));
z_mean = 1.0/n*@sum(row(i):z(i));
s = @sqrt(1.0/n*@sum(row(i):(z(i)-z_mean)^2));
@for(same(i,k):r(i,k)=@abs(z(i)-z(k)));
Rv = @max(same(i,k):r(i,k))+3;
d = @sum(same(i,k):(Rv-r(i,k))*@if((Rv-r(i,k))#ge#0,1,0));
@sum(col(j):(a(j))^2) = 1;
end

求解得到求解报告如下：