中國2012年數學建模大賽A題（葡萄與葡萄酒）+補充了pls迴歸係數的算法

grape=read.csv('d:/grape.csv') #葡萄數據
grape
grape=as.matrix(grape[,-1])
claret=read.csv('d:/claret.csv') #葡萄酒數據
claret
claret=as.matrix(claret[,-1])
x=grape
y=claret
x0=scale(x)
x0
y0=scale(y)
y0
m=t(x0)%*%y0%*%t(y0)%*%x0
m
eigen(m)
w1=eigen(m)$vectors[,1]
w1
v1=t(y0)%*%x0%*%w1/sqrt(as.matrix(eigen(m)$values)[1,])
v1
t1=x0%*%w1
t1 # 以上爲第一步(1)分別提取兩變量組的第一對成分，並使之相關性達最大。
u1=y0%*%v1
u1
library("pracma")
α1=inv(t(t1)%*%t1)%*%t(t1)%*%x0 #也可由t(x0)%*%t1/norm(t1,'2')^2算得α1(主成分t1與x0進行最小二乘迴歸求得的t1對每個x0的迴歸係數)
β1=inv(t(t1)%*%t1)%*%t(t1)%*%y0 #也可由t(y0)%*%t1/norm(t1,'2')^2算得β1(主成分t1與y0進行最小二乘迴歸求得的t1對每個y0的迴歸係數)
α1
t(α1)
β1
t(β1)
t(x0)%*%t1/norm(t1,'2')^2 # norm(t1,'2')爲svd(t1)即t1的最大奇異值,也可用sqrt(t(t1)%*%t1)求得
t(y0)%*%t1/norm(t1,'2')^2 # 以上爲第二步(2)創建y對T1的迴歸及x對T1的迴歸。
B=t(x0)%*%u1%*%inv(t(t1)%*%x0%*%t(x0)%*%u1)%*%t(t1)%*%y0 #保留第一對潛變量對應的標準化自變量x和標準化應變量y的pls迴歸係數矩陣(該矩陣公式參見‘KernelPartialLeastSquaresRegressionin Reproducing Kernel Hilbert Space’p102)

B
library("pls")
pls1=plsr(y0~x0,ncomp=1,validation='LOO',jackknife=T)
coef(pls1) #上式中B的求解等價於R的pls包中保留一個主成分的結果
px0=t1%*%α1
E1=x0-px0
py0=t1%*%β1
F1=y0-py0 # 以上爲第三步(3) 用殘差陣E1和F1代替x0和y0重複以上步驟。
m2=t(E1)%*%F1%*%t(F1)%*%E1
eigen(m2)
w2=eigen(m2)$vectors[,1]
w2
v2=t(F1)%*%E1%*%w2/sqrt(as.matrix(eigen(m2)$values)[1,])
v2
t2=E1%*%w2
t2
u2=F1%*%v2
u2
α2=inv(t(t2)%*%t2)%*%t(t2)%*%E1 #也可由t(E1)%*%t2/norm(t2,'2')^2算得α2
β2=inv(t(t2)%*%t2)%*%t(t2)%*%F1 #也可由t(F1)%*%t2/norm(t2,'2')^2算得β2
α2
β2
pE1=t2%*%α2
E2=E1-pE1
pF1=t2%*%β2
F2=F1-pF1
library("pls")
pls1=plsr(y0~x0,ncomp=10,validation='LOO',jackknife=T)
coef(pls1)
summary(pls1,what='all')
pls1=plsr(y0~x0,ncomp=3,validation='LOO',jackknife=T) #最終選擇3個主成分，從新創建模型
coef(pls1)
jack.test(pls1) #即對coef(pls1)生成的係數進行假設檢驗
scores(pls1) #即求t1=x0%*%w1
loadings(pls1) #即求α1
plot(pls1)
validationplot(pls1) #validationplot()函數能夠畫出PLS模型在不一樣主成分數下對應的RMSEP(由留一交叉驗證法算得的均方預測偏差根)
predict(pls1) #即求py0=t1%*%β1
predplot(pls1)