%基于樣條變換的偏最小二乘回歸
function plss()
x=[-5,-4,-3,-2,-1,0,1,2,3,4,5]';
y=[-19.649,-9.1919,-3.5907,4.5837,0.89779,-4.2108,-0.89676,-14.643,3.4229,2.8969,1.9406]';
f=zeros(11,7);
for n=1:11
    for l=0:6
        for k=0:4
            if x(n)>=(-5+(l-3+k)*2.5)
                f(n,l+1)=f(n,l+1)+(1/(factorial(3)*2.5^3))*(-1)^k*(factorial(4)/(factorial(4-k)*factorial(k)))*(x(n)-(-5+(l-3+k)*2.5))^3;
            else f(n,l+1)=f(n,l+1);
            end
        end
    end
end
[AA,aa,TT]=pls1(f,y,1);
a0=AA(1);
aa=AA(2:8);
s=-5:0.1:5;
ft=zeros(101,7);
for n=1:101
    for l=0:6
        ft(n,l+1)=0;
        t(n)=0;
        for k=0:4
            if s(n)>=(-5+(l-3+k)*2.5)
                ft(n,l+1)=ft(n,l+1)+(1/(factorial(3)*2.5^3))*(-1)^k*(factorial(4)/(factorial(4-k)*factorial(k)))*(s(n)-(-5+(l-3+k)*2.5))^3;
            else ft(n,l+1)=ft(n,l+1);
            end
        end
        t(n)=t(n)+aa(l+1)*ft(n,l+1);
    end
    t(n)=a0+ft(n,:)*aa';
    d(n)=0.1*(s(n)-1.2)^3;
end
plot(x,y,'r*',s,d,'r-',s,t,'b-'); grid on;
box on;
legend('樣本點','真實曲線','擬合曲線',4);
function [AA, aa, TT] = pls1(x,y,step)
[n,m]= size(x);
[yn,ym] = size(y);
aver = mean(x);
stdcov = std(x); 
E0 = ( x - aver(ones(n,1),:) )./ stdcov( ones(n,1),:);   
if ym == 1
    y1 = y;
else
    y1 = y(:,1);
end
F0 = (y1 - mean(y1) )./std(y1); 
E = E0;
F = F0;
for h = 1:step
    w = ( E'* F )/ norm( E'*F );
    W(:, h ) = w;
    t = E * w;
    T(:, h ) = t;  
    p = E'* t / ( norm( t)^2 );
    P(:, h ) = p;
    r = F'* t / ( norm( t)^2 );
    R( h ) = r;    
    E = E - t * p';
    F = F - t * r;
    if step == m & nargin == 2 
        if h == step     
            break;
        end
        SSi(h) = sum( (F0 - T*R').^2 );
        PRESSi1(h) = GetPRESS(E0,F0,h+1);
        if h>1
            if( 1- PRESSi1(h)/SSi(h-1)) < 0.0975
                break;
            end
        end
 end
end
TT = T;           
for i = 1: h   
    W_h_s = eye(m);
    for j = 1:(i-1)
        W_h_s = W_h_s * ( eye(m) - W(j) *( P(j) )');
    end
    W_h_s = W_h_s * W(:,i);  
    W_s(:,i) = W_h_s;
end
a_s =  R * W_s'; 
aa = a_s;              
ai = a_s * std(y1) ./ stdcov;
a0 = mean(y1) - sum( ( a_s * std(y1) ./ stdcov ) .* mean(x) );
AA = [a0,ai];
function re = GetPRESS(E0,F0,h)
[n,m] = size(E0);
F0n = size(F0);
press = 0;
for i = 1:n
    E = [ E0(1:i-1,:);  E0(i+1:n,:) ];
    F = [ F0(1:i-1,:);  F0(i+1:n,:) ];
 for j = 1:h         
        w = ( E'* F )/ norm( E'*F );
        W(:, j ) = w;
        t = E * w;
        T(:, j ) = t;    
        p = E'* t / ( norm( t)^2 );
        P(:, j ) = p;
        r = F'* t / ( norm( t)^2 );
        R( j ) = r;       
        E = E - t * p';
        F = F - t * r;
 end
 for k = 1: h
        W_h_s = eye(m);
        for j = 1:(k-1)
            W_h_s = W_h_s * ( eye(m) - W(j) *( P(j) )');
        end
        W_h_s = W_h_s * W(:,k);     
        W_s(:,k) = W_h_s;
 end
    F0i_s = E0(i,:) * (W_s * R');     
    F0i = F0( i);
    press = press + ( F0i - F0i_s )^2; 
end
re = press;