连续投影算法SPA的python实现

## 原理 连续投影算法（successive projections algorithm， SPA） 是前向特征变量选择方法。SPA利用向量的投影分析，通过将波长投影到其他波长上，比较投影向量大小，以投影向量最大的波长为待选波长，然后基于矫正模型选择最终的特征波长。SPA选择的是含有最少冗余信息及最小共线性的变量组合。该算法简要步骤如下。 记初始迭代向量为 $x_{k(0)}$，需要提取的变量个数为$N$,光谱矩阵为$J$列。 1. 任选光谱矩阵的1列（第$j$列），把建模集的第$j$列赋值给$x_j$，记为 $x_{k(0)}$。 2. 将未选入的列向量位置的集合记为$s$, $$ s=\lbrace j,1\leq{j}\leq{J}, j\notin \lbrace k(0), \cdots, k(n-1) \rbrace \rbrace $$ 3. 分别计算$x_j$对剩余列向量的投影： $$ P_{x_j} = x_j-(x^T_j x_{k(n-1)})x_{k(n-1)}(x^T_{k(n-1)}x_{k(n-1)})^{-1},j\in s $$ 4. 提取最大投影向量的光谱波长， $$ k(n) = arg(max(\| P_(x_j) \|), j \in s) $$ 5. 令$x_j = p_x, j \in s$。 6. $n = n + 1$，如果$n < N$,则按公式（1）循环计算。 最后，提取出的变量为$\lbrace x_{k(n)} = 0, \cdots, N-1 \rbrace$。对应每一次循环中的$k(0)$和$N$，分别建立多元线性回归分析（MLR）模型，得到建模集交互验证均方根误差（RMSECV），对应不同的候选子集，其中最小的RMSECV值对应的$k(0)$和$N$就是最优值。一般SPA选择的特征波长分数$N$不能很大。 ​ -------------------------摘自《光谱及成像技术在农业中的应用》P130 ## 快速使用 ## 0. WARNING 当时python开发环境主要版本 python == 3.7 Scikit-learn == 0.24.0 后续版本可能报错 建议使用虚拟环境运行该脚本 脚本按照MATLAB版本一比一复刻 仅供参考 可按实际情况进行修改 如果您编程能力较弱 对python理解不足 不建议参考使用本项目代码 本项目代码未经大量测试 无法保证通用性。 如果您对本项目进行了结构性更改，请一定通知我进行同步更改。 ### 1. 数据读入 ```python # 导入 pandas 读取数据 import pandas as pd import numpy as np # 读取数据 data = pd.read_csv("./data/peach_spectra_brix.csv") data[:5] ``` ### 2. 数据分离 ```python # m * n print("数据矩阵 data.shape：",data.shape) # 50个样本， 600个 波段 第一列是 桃子糖度值 需要分离开 X = data.values[:,1:] # 等同操作 #X = data.drop(['Brix'], axis=1) y = data.values[:,0] # 等同操作 # y = data.loc[:,'Brix'].values print(f"X.shape:{X.shape}, y.shape:{y.shape}") X,y ``` ### 3. 导入SPA ```python # 导入SPA包 import SPA # 导入spa对象 spa = SPA.SPA() # 数据归一化 from sklearn.preprocessing import MinMaxScaler min_max_scaler = MinMaxScaler(feature_range=(-1, 1)) # 这里feature_range根据需要自行设置，默认（0,1） X_ = min_max_scaler.fit_transform(X) print(X[1,:5]) print(X_[1,:5]) # 建模集测试集分割 from sklearn.model_selection import train_test_split # 注意 X_ 分割 # 若存在 运行后出现 波段选择为 最小值 可适当调整 建模集与测试集比例 test_size 值 0.3 - 0.5 Xcal, Xval, ycal, yval = train_test_split(X_, y, test_size=0.4, random_state=0) Xcal.shape,Xval.shape ``` ### 4. 建模筛选 ```python # 建模筛选 # m_max 默认为 50(Xcal样本大于52) ,如果 Xcal(m*n) m < 50 m_max=m-2 var_sel, var_sel_phase2 = spa.spa( Xcal, ycal, m_min=2, m_max=28,Xval=Xval, yval=yval, autoscaling=1) ``` ### 5. 导出波段 ```python # 导出 筛选光谱波段 # spa 返回的是 列号 并不是光谱数据 # 获取 波段列表 absorbances = data.columns.values[1:] print("波段(前5个)",absorbances[:5]) # spa 筛选出的波段 以 test_size=0.3 为例 print("spa 筛选出的波段(以 test_size=0.3 为例):",absorbances[var_sel]) # 导出 筛选波段光谱数据 X_select = X[:,var_sel] print("X_select.shape:",X_select.shape) ``` [示例&建模参考](https://gitee.com/aBugsLife/SPA/blob/master/SPA%E7%A4%BA%E4%BE%8B.ipynb) ## 注意事项 1. 光谱矩阵（m * n）**m行为样本，n列为波段** 2. 进行建模前需要对光谱进行 建模集测试集分割 与 数据归一化 ，可先进行分割再归一，也可以先归一再分割，下边为分割再归一 ````python from sklearn.model_selection import train_test_split from sklearn.preprocessing import MinMaxScaler Xcal, Xval, ycal, yval = train_test_split(x, y, test_size=0.4, random_state=0) min_max_scaler = MinMaxScaler(feature_range=(-1, 1)) Xcal = min_max_scaler.fit_transform(Xcal) Xval = min_max_scaler.transform(Xval) ```` 示例数据来源：[nirpyresearch.com](https://nirpyresearch.com/)