get_derivative.m

来自「主成分分析和偏最小二乘SquaresPrincipal成分分析（ PCA ）和偏」· M 代码 · 共 107 行

function [Yr,state]=get_derivative(method,order,X,Y,Xr);
% method - now available only 'poly^4'
% order- order of derivative
% given: Y=f(X)
% for order=1 function compute Yr=DY/DX=f'(Xr)
% for order=2 function compute Yr=D2Y/DX2=f''(Xr)
% for order=0 function compute Yr=f(Xr) using poly^4 interpolation
% state-{1-OK, 0-Error}

% last modified 4.10.04

state=0;
switch method
    case 'poly^4' 
        Yr=ones(size(Xr));
        LR=length(Xr);
        MaxL=length(X);
        MinL=1;
        for i=1:LR
   %        R=find(X(MinL:MaxL)<=Xr(i));
            tp=MinL;
            R=0;
            while tp<=MaxL
                if X(tp)<=Xr(i)
                    R=R+1;
                else
                    break
                end
                tp=tp+1;
            end
                
            if R==0
                return
            end
            Cp=max(R)+MinL-1;
            if Cp<3
                if Cp==2
                    XJ=X(Cp-1:Cp+3);
                    YJ=Y(Cp-1:Cp+3);
                else
                    XJ=X(Cp:Cp+4);
                    YJ=Y(Cp:Cp+4);
                end
            elseif Cp>MaxL-2
                if Cp==MaxL-1
                    XJ=X(Cp-3:Cp+1);
                    YJ=Y(Cp-3:Cp+1);
                else
                    XJ=X(Cp-4:Cp);
                    YJ=Y(Cp-4:Cp);
                end
            else
                XJ=X(Cp-2:Cp+2);
                YJ=Y(Cp-2:Cp+2);
            end
              for deg=4:-1:2
%                  deg=4;
                [p,ok] = lab432_polyfit(XJ,YJ,deg);
                  if ok
                      break
                  end
            end
            switch order
                case 0
                    Yr(i)=polyval(p,Xr(i));
                case 1
                    Yr(i)=polyval(polyder(p),Xr(i));
                case 2
                    Yr(i)=polyval(polyder(polyder(p)),Xr(i));
                otherwise
                    disp('Invalid order. Order must be 0,1 or 2');
            end
            MinL=Cp;
        end
end
state=1;               


function [p,ok] = lab432_polyfit(x,y,n)
%POLYFIT Fit polynomial to data.
%   POLYFIT(X,Y,N) finds the coefficients of a polynomial P(X) of
%   degree N that fits the data, P(X(I))~=Y(I), in a least-squares sense.

ok=1;
x = x(:);
y = y(:);
if ~isequal(size(x),size(y))
    error('X and Y vectors must be the same size.')
end


V(:,n+1) = ones(length(x),1);
for j = n:-1:1
   V(:,j) = x.*V(:,j+1);
end

% Solve least squares problem
[Q,R] = qr(V,0);
ws = warning('off','all'); 
p = R\(Q'*y);    % Same as p = V\y;
warning(ws);
if size(R,2) > size(R,1)
   ok=0;
elseif condest(R) > 1.0e10
    ok=0;
end
p = p.';