第9讲--多元线性回归--主成分回归

程序代码—预报
对新测定Snew: C=np.dot(Snew, Alast)
扩展--能否用MLR、PCA类PCR类
• 传递X，Y给PCR
• PCR内，以X调用PCA，确定主成分数
• 根据确定的主成分数，确定T、P，以T，Y建模，并结合
P确定回归系数 • 建立预报方法。
扩展--能否用MLR、PCA类PCR类
T和P都是列正交矩阵
T的第i列ti的模，就是第i个特征值λi E为残差矩阵，对应噪声
11
主成分示例
方差最大方向
NIPALS算法每次只求一个主成分，目前最大散差方向
12
仪器的信噪比
仪器测量时，信号强度要远远大于噪声
信号的数据的方差要远远大于噪声的方差 所以，PCA可以区别噪声
第9讲 多元线性回归 主成分及其回归
多元线性回归解决的问题
系数矩阵 Y=XA 建模：求解回归系数A，该过程称为建
模 预报：在A已知时，对于新测Xnew， 预报Ynew，称为预报
例子
某保健品含片产品，说明书标明：由营养
物质A、B、C组成，产品标注中写出了每 片中A、B、C物质的含量。问，如何认定 ？
矩阵中有3个有效特征值 根据有效特征值，设定PCA 的得分和载荷
实例—光谱矩阵的SVD分解
根据主成分，规划得分U和载荷矩阵P SVD ：X=USV PCA：X=TPt T=US ，P=Vt i=len(lamda) S=np.zeros ((i,i)) S[:i,:i]=np.diag (lamda) T = np.dot (U,S) V=V.T P =V T = T[:,:k] P = P[:,:k]
程序代码—建模过程
P=V.T for i in range(len(lmada)-1): temp = lamda[i]/lamda[i+1] print (temp) k=int(input(“主成分数为:”)) T = T[:,:k] P = P[:,:k] TtT=np.dot(T.T,T) inv = np.linalg.inv(TtT) A=np.dot(inv, T.T) Alast=np.dot(P,A)
主成分回归
数据 E:\学校教学\python\S-093843.txt E:\学校教学\python\C-093843.txt 求解方程 C=SA S是6*16的矩阵 所以StS的逆不存在 S 是光谱矩阵，光谱的不同波长间线性相关，所以 可以用PCR
程序代码—建模过程
S=np.mafromtxt(“E:\\学校教学\\python\\S-093843.txt ") S=S.data
通过特征值比值判断有效变量数
在λi/ λi+i，应该达到最大值 根据i值，取T和P的前i列，即可扔掉噪声
15
主成分回归PCR
Principle Component Regression 是多元线性回归！ 原来 Y=XA 现在 Y=TA T为X的主成分得分，即X经PCA分解后的得分 因为 T只是 X 的线性组合，提取了线性相关的 部分，且只取前i列，所以模型稳定，去掉噪 声
配置A、B、C的一组溶液，建立浓度与光
吸收的关系。既建模求回归系数 将药片配置成溶液，测吸光，利用上面的 模型，预报浓度。
建模公式推导
Y=XA XtY=XtXA (XtX)-1XtY=A
E:\学校教学\python\X.txt E:\学校教学\python\Y.txt
问题求解的关键步骤是什么？
样例
原 数 据
x y
0.9 1.1 1.2 1.0
0.8
0.87 2
2.2
1.9 2.1 1.7 2.5
0.92 1.1
1.81 1.9
图
PCA后
t1 t2 1.486 -0.208 1.485 0.075 1.216 -0.081 1.393 -0.158 2.694 0.142 2.898 0.221 2.545 0.149 3.253 -0.273
print("相邻特征值比值") compare = pcr.confirmPCs() print(compare) k=int(input("确定主成分数：")) pcr.model(k) ans=pcr.predict(S)
方程数与未知数的关系
设有规律上符合如下方程的一 组实验数据
y= ax+b
通过实验，不断变更x，测得对应的y 求a,b的值，需要几组这样的数据？ 唯一解 最小二乘解
y1 y2 … yn
=
x1 x2 … பைடு நூலகம்n
1 1 1 1
矩阵形式
a b XtX是2*2的矩阵
方程数与未知数的关系
设有规律上符合如下方程的一 组实验数据
16
numpy中主成分分解—SVD分解
实矩阵的SVD(Singular Value Decomposition，奇
异值分解 )分解： 分解结果：A=USV 其中S是对角矩阵
numpy中主成分分解---SVD
程序代码： B = np.linalg.svd(A,full_matrices=False) full_matrices=False一定要写，否则会按复数分解 分解结果： U=B[0] lamda=B[1] V = B[2] Lamda是所有的特征值，可以计算相邻比值，决定 主成分，它不是一个矩阵
y= 1.2 x1 + 0.9 x2 + 3.3 x3
通过实验，不断变更x1、x2、x3，测得对应的y
需要几组这样的数据？
唯一解
最小二乘解
方程数小于未知数，一定无解吗
如果x个数很多，样本打不到要求，怎么办 ？ y= 1.2 x1 + 0.9 x2 + 3.3 x3
当X1,X2,X3存在线性相关时，问题会怎样？
PCA类
self.T = np.dot (U,S) V=V.T self.P = V compare=[] for i in range(len(lamda)-1): temp = lamda[i]/lamda[i+1] compare.append(temp) return U,S,V,compare def PCAdecompose(self,k): T = self.T[:,:k] P = self.P[:,:k] return T,P
实例—光谱矩阵的SVD分解
数据：E:\学校教学\教改项目教材\数据\S-093790.txt 是一个16*6的矩阵 看看能求解个特征值？16个？ 6个？96个？
实例—光谱矩阵的SVD分解
data=np.mafromtxt("E:\\学校教学\\教改项目教材\\ 数据\\S-093790.txt") data=data.data B = np.linalg.svd(data,full_matrices=False) >>> B[1] array([ 5.48250094e+00, 1.10440342e+00, 3.27012276e-01, 3.23153080e-03, 2.19720845e-03, 1.11546885e03])
可否编写PCA类
传递矩阵给类 求得T、P矩阵，特征值比值列表 根据特征值比值，规划T和P
PCA类
import numpy as np class PCA: def __init__(self, A): self.A=A def SVDdecompose(self): B = np.linalg.svd(self.A,full_matrices=False) U=B[0] lamda=B[1] V = B[2] i=len(lamda) S=np.zeros ((i,i)) S[:i,:i]=np.diag (lamda)
import numpy as np from PCA import PCA from MLR import MLR class PCR: def __init__(self,X,Y): self.X=X self.Y=Y def confirmPCs(self): pca=PCA(self.X) U,S,V,compare=pca.SVDdecompose() return compare
C=np.mafromtxt(“E:\\学校教学\\python\\C-093843.txt ") C=C.data B = np.linalg.svd(data,full_matrices=False) U=B[0] lamda=B[1] i=len(lamda) S=np.zeros ((i,i)) S[:i,:i]=np.diag (lamda) T = np.dot (U,S) V=B[2]
调用函数
S=np.mafromtxt("E:\\学校教学\\python\\S-093843.txt") S=S.data
S=S.T 要一行一个样本，所以转置 C=np.mafromtxt("E:\\学校教学\\python\\C-093843.txt") C=C.data C=C.T pcr=PCR(S,C)
PCA类调用
应该先调用decompose方法，根据返回的
特征之比值，确定主成分 再调用PCAdecompose方法，设定得分和 载荷矩阵
主成分回归
原来 Y=XA 现在 Y=TA
TtY=TtTA TtT-1TtY=A 求得最终的回归系数：主成分 X=TPt 因为P是正交矩阵，所以 T=XP Y=TA=XPA=XAnew
扩展----能否用MLR、PCA类
def model(self,PCs): pca=PCA(self.X) U,S,V,compare=pca.SVDdecompose() T,P=pca.PCAdecompose(PCs) mlr=MLR(T,self.Y) mlr.modelling() self.A=np.dot(P,mlr.A) def predict(self,Xnew): ans = np.dot(Xnew,self.A) return ans