精华内容
下载资源
问答
  • 偏最小二乘法Matlab源代码,自己编写
  • 一个偏最小二乘法的应用实例,希望能对读者有帮助
  • matlab实现偏最小二乘法。。。。。 %% Principal Component Analysis and Partial Least Squares % Principal Component Analysis (PCA) and Partial Least Squares (PLS) are % widely used tools. This code is to...
  • 偏最小二乘法 matlab程序

    千次阅读 2010-01-07 13:22:00
    偏最小二乘法 matlab程序函数:function [T,P,U,Q,B,W] = pls (X,Y,tol2)% PLS Partial Least Squares Regrassion%% [T,P,U,Q,B,Q] = pls(X,Y,tol) performs particial least squares regrassion% between the ...

    偏最小二乘法 matlab程序

    函数:

    function [T,P,U,Q,B,W] = pls (X,Y,tol2)
    % PLS   Partial Least Squares Regrassion
    %
    % [T,P,U,Q,B,Q] = pls(X,Y,tol) performs particial least squares regrassion
    % between the independent variables, X and dependent Y as
    % X = T*P' + E;
    % Y = U*Q' + F = T*B*Q' + F1;
    %
    % Inputs:
    % X     data matrix of independent variables
    % Y     data matrix of dependent variables
    % tol   the tolerant of convergence (defaut 1e-10)
    %
    % Outputs:
    % T     score matrix of X
    % P     loading matrix of X
    % U     score matrix of Y
    % Q     loading matrix of Y
    % B     matrix of regression coefficient
    % W     weight matrix of X
    %
    % Using the PLS model, for new X1, Y1 can be predicted as
    % Y1 = (X1*P)*B*Q' = X1*(P*B*Q')
    % or
    % Y1 = X1*(W*inv(P'*W)*inv(T'*T)*T'*Y)
    %
    % Without Y provided, the function will return the principal components as
    % X = T*P' + E
    %
    % Example: taken from Geladi, P. and Kowalski, B.R., "An example of 2-block
    % predictive partial least-squares regression with simulated data",
    % Analytica Chemica Acta, 185(1996) 19--32.
    %{
    x=[4 9 6 7 7 8 3 2;6 15 10 15 17 22 9 4;8 21 14 23 27 36 15 6;
    10 21 14 13 11 10 3 4; 12 27 18 21 21 24 9 6; 14 33 22 29 31 38 15 8;
    16 33 22 19 15 12 3 6; 18 39 26 27 25 26 9 8;20 45 30 35 35 40 15 10];
    y=[1 1;3 1;5 1;1 3;3 3;5 3;1 5;3 5;5 5];
    % leave the last sample for test
    N=size(x,1);
    x1=x(1:N-1,:);
    y1=y(1:N-1,:);
    x2=x(N,:);
    y2=y(N,:);
    % normalization
    xmean=mean(x1);
    xstd=std(x1);
    ymean=mean(y1);
    ystd=std(y);
    X=(x1-xmean(ones(N-1,1),:))./xstd(ones(N-1,1),:);
    Y=(y1-ymean(ones(N-1,1),:))./ystd(ones(N-1,1),:);
    % PLS model
    [T,P,U,Q,B,W]=pls(X,Y);
    % Prediction and error
    yp = (x2-xmean)./xstd * (P*B*Q');
    fprintf('Prediction error: %g/n',norm(yp-(y2-ymean)./ystd));
    %}
    %
    % By Yi Cao at Cranfield University on 2nd Febuary 2008
    %
    % Reference:
    % Geladi, P and Kowalski, B.R., "Partial Least-Squares Regression: A
    % Tutorial", Analytica Chimica Acta, 185 (1986) 1--7.
    %

    % Input check
    error(nargchk(1,3,nargin));
    error(nargoutchk(0,6,nargout));
    if nargin<2
        Y=X;
    end
    tol = 1e-10;
    if nargin<3
        tol2=1e-10;
    end

    % Size of x and y
    [rX,cX]  =  size(X);
    [rY,cY]  =  size(Y);
    assert(rX==rY,'Sizes of X and Y mismatch.');

    % Allocate memory to the maximum size
    n=max(cX,cY);
    T=zeros(rX,n);
    P=zeros(cX,n);
    U=zeros(rY,n);
    Q=zeros(cY,n);
    B=zeros(n,n);
    W=P;
    k=0;
    % iteration loop if residual is larger than specfied
    while norm(Y)>tol2 && k<n
        % choose the column of x has the largest square of sum as t.
        % choose the column of y has the largest square of sum as u.   
        [dummy,tidx] =  max(sum(X.*X));
        [dummy,uidx] =  max(sum(Y.*Y));
        t1 = X(:,tidx);
        u = Y(:,uidx);
        t = zeros(rX,1);

        % iteration for outer modeling until convergence
        while norm(t1-t) > tol
            w = X'*u;
            w = w/norm(w);
            t = t1;
            t1 = X*w;
            q = Y'*t1;
            q = q/norm(q);
            u = Y*q;
        end
        % update p based on t
        t=t1;
        p=X'*t/(t'*t);
        pnorm=norm(p);
        p=p/pnorm;
        t=t*pnorm;
        w=w*pnorm;
       
        % regression and residuals
        b = u'*t/(t'*t);
        X = X - t*p';
        Y = Y - b*t*q';
       
        % save iteration results to outputs:
        k=k+1;
        T(:,k)=t;
        P(:,k)=p;
        U(:,k)=u;
        Q(:,k)=q;
        W(:,k)=w;
        B(k,k)=b;
        % uncomment the following line if you wish to see the convergence
    %     disp(norm(Y))
    end
    T(:,k+1:end)=[];
    P(:,k+1:end)=[];
    U(:,k+1:end)=[];
    Q(:,k+1:end)=[];
    W(:,k+1:end)=[];
    B=B(1:k,1:k);

     

     

    例子:——————————————————————————————————————————————————————————

    %% Principal Component Analysis and Partial Least Squares
    % Principal Component Analysis (PCA) and Partial Least Squares (PLS) are
    % widely used tools. This code is to show their relationship through the
    % Nonlinear Iterative PArtial Least Squares (NIPALS) algorithm.

    %% The Eigenvalue and Power Method
    % The NIPALS algorithm can be derived from the Power method to solve the
    % eigenvalue problem. Let x be the eigenvector of a square matrix, A,
    % corresponding to the eignvalue s:
    %
    % $$Ax=sx$$
    %
    % Modifying both sides by A iteratively leads to
    %
    % $$A^kx=s^kx$$
    %
    % Now, consider another vectro y, which can be represented as a linear
    % combination of all eigenvectors:
    %
    % $$y=/sum_i^n b_ix_i=Xb$$
    %
    % where
    %
    % $$X=/left[x_1/,/,/, /cdots/,/,/, x_n /right]$$
    %
    % and
    %
    % $$b = /left[b_1/,/,/, /cdots/,/,/, b_n /right]^T$$
    %
    % Modifying y by A gives
    %
    % $$Ay=AXb=XSb$$
    %
    % Where S is a diagnal matrix consisting all eigenvalues. Therefore, for
    % a large enough k,
    %
    % $$A^ky=XS^kb/approx /alpha x_1$$
    %
    % That is the iteration will converge to the direction of x_1, which is the
    % eigenvector corresponding to the eigenvalue with the maximum module.
    % This leads to the following Power method to solve the eigenvalue problem.

    A=randn(10,5);
    % sysmetric matrix to ensure real eigenvalues
    B=A'*A;
    %find the column which has the maximum norm
    [dum,idx]=max(sum(A.*A));
    x=A(:,idx);
    %storage to judge convergence
    x0=x-x;
    %convergence tolerant
    tol=1e-6;
    %iteration if not converged
    while norm(x-x0)>tol
        %iteration to approach the eigenvector direction
        y=A'*x;
        %normalize the vector
        y=y/norm(y);
        %save previous x
        x0=x;
        %x is a product of eigenvalue and eigenvector
        x=A*y;
    end
    % the largest eigen value corresponding eigenvector is y
    s=x'*x;
    % compare it with those obtained with eig
    [V,D]=eig(B);
    [d,idx]=max(diag(D));
    v=V(:,idx);
    disp(d-s)
    % v and y may be different in signs
    disp(min(norm(v-y),norm(v+y)))

    %% The NIPALS Algorithm for PCA
    % The PCA is a dimension reduction technique, which is based on the
    % following decomposition:
    %
    % $$X=TP^T+E$$
    %
    % Where X is the data matrix (m x n) to be analysed, T is the so called
    % score matrix (m x a), P the loading matrix (n x a) and E the residual.
    % For a given tolerance of residual, the number of principal components, a,
    % can be much smaller than the orginal variable dimension, n.
    % The above power algorithm can be extended to get T and P by iteratively
    % subtracting A (in this case, X) by x*y' (in this case, t*p') until the
    % given tolerance satisfied. This is the so called NIPALS algorithm.

    % The data matrix with normalization
    A=randn(10,5);
    meanx=mean(A);
    stdx=std(A);
    X=(A-meanx(ones(10,1),:))./stdx(ones(10,1),:);
    B=X'*X;
    % allocate T and P
    T=zeros(10,5);
    P=zeros(5);
    % tol for convergence
    tol=1e-6;
    % tol for PC of 95 percent
    tol2=(1-0.95)*5*(10-1);
    for k=1:5
        %find the column which has the maximum norm
        [dum,idx]=max(sum(X.*X));
        t=A(:,idx);
        %storage to judge convergence
        t0=t-t;
        %iteration if not converged
        while norm(t-t0)>tol
            %iteration to approach the eigenvector direction
            p=X'*t;
            %normalize the vector
            p=p/norm(p);
            %save previous t
            t0=t;
            %t is a product of eigenvalue and eigenvector
            t=X*p;
        end
        %subtracing PC identified
        X=X-t*p';
        T(:,k)=t;
        P(:,k)=p;
        if norm(X)<tol2
            break
        end
    end
    T(:,k+1:5)=[];
    P(:,k+1:5)=[];
    S=diag(T'*T);
    % compare it with those obtained with eig
    [V,D]=eig(B);
    [D,idx]=sort(diag(D),'descend');
    D=D(1:k);
    V=V(:,idx(1:k));
    fprintf('The number of PC: %i/n',k);
    fprintf('norm of score difference between EIG and NIPALS: %g/n',norm(D-S));
    fprintf('norm of loading difference between EIG and NIPALS: %g/n',norm(abs(V)-abs(P)));

    %% The NIPALS Algorithm for PLS
    % For PLS, we will have two sets of data: the independent X and dependent
    % Y. The NIPALS algorithm can be used to decomposes both X and Y so that
    %
    % $$X=TP^T+E,/,/,/,/,Y=UQ^T+F,/,/,/,/,U=TB$$
    %
    % The regression, U=TB is solved through least sequares whilst the
    % decompsition may not include all components. That is why the approach is
    % called partial least squares. This algorithm is implemented in the PLS
    % function.

    %% Example: Discriminant PLS using the NIPALS Algorithm
    % From Chiang, Y.Q., Zhuang, Y.M and Yang, J.Y, "Optimal Fisher
    % discriminant analysis using the rank decomposition", Pattern Recognition,
    % 25 (1992), 101--111.

    % Three classes data, each has 50 samples and 4 variables.
    x1=[5.1 3.5 1.4 0.2; 4.9 3.0 1.4 0.2; 4.7 3.2 1.3 0.2; 4.6 3.1 1.5 0.2;...
        5.0 3.6 1.4 0.2; 5.4 3.9 1.7 0.4; 4.6 3.4 1.4 0.3; 5.0 3.4 1.5 0.2; ...
        4.4 2.9 1.4 0.2; 4.9 3.1 1.5 0.1; 5.4 3.7 1.5 0.2; 4.8 3.4 1.6 0.2; ...
        4.8 3.0 1.4 0.1; 4.3 3.0 1.1 0.1; 5.8 4.0 1.2 0.2; 5.7 4.4 1.5 0.4; ...
        5.4 3.9 1.3 0.4; 5.1 3.5 1.4 0.3; 5.7 3.8 1.7 0.3; 5.1 3.8 1.5 0.3; ...
        5.4 3.4 1.7 0.2; 5.1 3.7 1.5 0.4; 4.6 3.6 1.0 0.2; 5.1 3.3 1.7 0.5; ...
        4.8 3.4 1.9 0.2; 5.0 3.0 1.6 0.2; 5.0 3.4 1.6 0.4; 5.2 3.5 1.5 0.2; ...
        5.2 3.4 1.4 0.2; 4.7 3.2 1.6 0.2; 4.8 3.1 1.6 0.2; 5.4 3.4 1.5 0.4; ...
        5.2 4.1 1.5 0.1; 5.5 4.2 1.4 0.2; 4.9 3.1 1.5 0.2; 5.0 3.2 1.2 0.2; ...
        5.5 3.5 1.3 0.2; 4.9 3.6 1.4 0.1; 4.4 3.0 1.3 0.2; 5.1 3.4 1.5 0.2; ...
        5.0 3.5 1.3 0.3; 4.5 2.3 1.3 0.3; 4.4 3.2 1.3 0.2; 5.0 3.5 1.6 0.6; ...
        5.1 3.8 1.9 0.4; 4.8 3.0 1.4 0.3; 5.1 3.8 1.6 0.2; 4.6 3.2 1.4 0.2; ...
        5.3 3.7 1.5 0.2; 5.0 3.3 1.4 0.2];
    x2=[7.0 3.2 4.7 1.4; 6.4 3.2 4.5 1.5; 6.9 3.1 4.9 1.5; 5.5 2.3 4.0 1.3; ...
        6.5 2.8 4.6 1.5; 5.7 2.8 4.5 1.3; 6.3 3.3 4.7 1.6; 4.9 2.4 3.3 1.0; ...
        6.6 2.9 4.6 1.3; 5.2 2.7 3.9 1.4; 5.0 2.0 3.5 1.0; 5.9 3.0 4.2 1.5; ...
        6.0 2.2 4.0 1.0; 6.1 2.9 4.7 1.4; 5.6 2.9 3.9 1.3; 6.7 3.1 4.4 1.4; ...
        5.6 3.0 4.5 1.5; 5.8 2.7 4.1 1.0; 6.2 2.2 4.5 1.5; 5.6 2.5 3.9 1.1; ...
        5.9 3.2 4.8 1.8; 6.1 2.8 4.0 1.3; 6.3 2.5 4.9 1.5; 6.1 2.8 4.7 1.2; ...
        6.4 2.9 4.3 1.3; 6.6 3.0 4.4 1.4; 6.8 2.8 4.8 1.4; 6.7 3.0 5.0 1.7; ...
        6.0 2.9 4.5 1.5; 5.7 2.6 3.5 1.0; 5.5 2.4 3.8 1.1; 5.5 2.4 3.7 1.0; ...
        5.8 2.7 3.9 1.2; 6.0 2.7 5.1 1.6; 5.4 3.0 4.5 1.5; 6.0 3.4 4.5 1.6; ...
        6.7 3.1 4.7 1.5; 6.3 2.3 4.4 1.3; 5.6 3.0 4.1 1.3; 5.5 2.5 5.0 1.3; ...
        5.5 2.6 4.4 1.2; 6.1 3.0 4.6 1.4; 5.8 2.6 4.0 1.2; 5.0 2.3 3.3 1.0; ...
        5.6 2.7 4.2 1.3; 5.7 3.0 4.2 1.2; 5.7 2.9 4.2 1.3; 6.2 2.9 4.3 1.3; ...
        5.1 2.5 3.0 1.1; 5.7 2.8 4.1 1.3];
    x3=[6.3 3.3 6.0 2.5; 5.8 2.7 5.1 1.9; 7.1 3.0 5.9 2.1; 6.3 2.9 5.6 1.8; ...
        6.5 3.0 5.8 2.2; 7.6 3.0 6.6 2.1; 4.9 2.5 4.5 1.7; 7.3 2.9 6.3 1.8; ...
        6.7 2.5 5.8 1.8; 7.2 3.6 6.1 2.5; 6.5 3.2 5.1 2.0; 6.4 2.7 5.3 1.9; ...
        6.8 3.0 5.5 2.1; 5.7 2.5 5.0 2.0; 5.8 2.8 5.1 2.4; 6.4 3.2 5.3 2.3; ...
        6.5 3.0 5.5 1.8; 7.7 3.8 6.7 2.2; 7.7 2.6 6.9 2.3; 6.0 2.2 5.0 1.5; ...
        6.9 3.2 5.7 2.3; 5.6 2.8 4.9 2.0; 7.7 2.8 6.7 2.0; 6.3 2.7 4.9 1.8; ...
        6.7 3.3 5.7 2.1; 7.2 3.2 6.0 1.8; 6.2 2.8 4.8 1.8; 6.1 3.0 4.9 1.8; ...
        6.4 2.8 5.6 2.1; 7.2 3.0 5.8 1.6; 7.4 2.8 6.1 1.9; 7.9 3.8 6.4 2.0; ...
        6.4 2.8 5.6 2.2; 6.3 2.8 5.1 1.5; 6.1 2.6 5.6 1.4; 7.7 3.0 6.1 2.3; ...
        6.3 3.4 5.6 2.4; 6.4 3.1 5.5 1.8; 6.0 3.0 4.8 1.8; 6.9 3.1 5.4 2.1; ...
        6.7 3.1 5.6 2.4; 6.9 3.1 5.1 2.3; 5.8 2.7 5.1 1.9; 6.8 3.2 5.9 2.3; ...
        6.7 3.3 5.7 2.5; 6.7 3.0 5.2 2.3; 6.3 2.5 5.0 1.9; 6.5 3.0 5.2 2.0; ...
        6.2 3.4 5.4 2.3; 5.9 3.0 5.1 1.8];
    %Split data set into training (1:25) and testing (26:50)
    idxTrain = 1:25;
    idxTest = 26:50;

    % Combine training data with normalization
    X = [x1(idxTrain,:);x2(idxTrain,:);x3(idxTrain,:)];
    % Define class indicator as Y
    Y = kron(eye(3),ones(25,1));
    % Normalization
    xmean = mean(X);
    xstd = std(X);
    ymean = mean(Y);
    ystd = std(Y);
    X = (X - xmean(ones(75,1),:))./xstd(ones(75,1),:);
    Y = (Y - ymean(ones(75,1),:))./ystd(ones(75,1),:);
    % Tolerance for 90 percent score
    tol = (1-0.9) * 25 * 4;
    % Perform PLS
    [T,P,U,Q,B] = pls(X,Y,tol);
    % Results
    fprintf('Number of components retained: %i/n',size(B,1))
    % Predicted classes
    X1 = (x1(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
    X2 = (x2(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
    X3 = (x3(idxTest,:) - xmean(ones(25,1),:))./xstd(ones(25,1),:);
    Y1 = X1 * (P*B*Q');
    Y2 = X2 * (P*B*Q');
    Y3 = X3 * (P*B*Q');
    Y1 = Y1 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
    Y2 = Y2 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
    Y3 = Y3 .* ystd(ones(25,1),:) + ymean(ones(25,1),:);
    % Class is determined from the cloumn which is most close to 1
    [dum,classid1]=min(abs(Y1-1),[],2);
    [dum,classid2]=min(abs(Y2-1),[],2);
    [dum,classid3]=min(abs(Y3-1),[],2);
    bar(1:25,classid1,'b');
    hold on
    bar(26:50,classid2,'r');
    bar(51:75,classid3,'g');
    hold off

    % The results show that most samples are classified correctly.

    展开全文
  • 详细的讲解非线性偏最小二乘法的运算过程,可以将其直接用来计算
  • 偏最小二乘法Matlab源码(2008-09-21 09:31:21) 标签 杂谈? ? 所谓偏最小二乘法就是指在做基于最小二乘法的线性回归分析之前对数据集进行主成分分析降维下面的源码是没有删减的GreenSim团队免费提供您使用转载请...
  • matlab计算偏最小二乘法的程序-matlab计算偏最小二乘法的程序.rar 偏最小二乘法 1建模原理 2计算方法推导 3交叉有效性 附录(源程序)
  • 偏最小二乘法代码PLS_connectome 此存储库包含对连接数据 (fMRI) 运行偏最小二乘 (PLS) 分析的脚本 用法 McIntosh 等人 (1996) 的 pls Matlab 工具箱中的函数必须在 Matlab 环境中可用。 当前脚本格式化功能连接矩阵...
  • 很好的matlab代码 偏最小二乘法
  • pls偏最小二乘法matlab实现,网上下载的,都打包到一起了,慢慢看吧。我从中找到了自己想要的,希望对你有帮助!
  • 全秩和部分秩偏最小二乘回归。 这些函数接受 X 和 Y 数据作为矩阵或表格。 X 和 Y 被转换为矩阵并进行处理。 PLS 组件的最佳数量是通过留一法交叉验证找到的,但它可以由用户修改。 输出被组织在用户定义的结构变量...
  • 偏最小二乘法PLS(matlab自带代码)

    千次阅读 2020-08-25 22:43:05
    偏最小二乘法(NIPALS经典实现--未简化) 偏最小二乘法 基本性质推导 偏最小二乘法(SIMPLS---未简化) 前面讲了许多PLS的性质,这里介绍一下PLS的matlab自带的PLS代码,对PLS有一个全面的认识 function ...
    1. 偏最小二乘法(NIPALS经典实现--未简化)
    2. 偏最小二乘法 基本性质推导
    3. 偏最小二乘法(SIMPLS---未简化)

     

    前面讲了许多PLS的性质,这里介绍一下PLS的matlab自带的PLS代码,对PLS有一个全面的认识

    matlab自带的plsregress函数实现了SIMPLS算法,主要是对X0和Y0的协方差矩阵不断做正交分解

    function [Xloadings,Yloadings,Xscores,Yscores, ...
                        beta,pctVar,mse,stats] = plsregress(X,Y,ncomp,varargin)
    %PLSREGRESS Partial least squares regression.
    
    
    ....
    % Center both predictors and response, and do PLS
    %对X,Y进行中心化,在模型中去掉截距项
    meanX = mean(X,1);
    meanY = mean(Y,1);
    X0 = bsxfun(@minus, X, meanX);
    Y0 = bsxfun(@minus, Y, meanY);
    
    % 根据输出的需要,返回不同的参数,下面重点关注SIMPLS,
    %这是Dejone在九几年提出的,matlab用这个算法来实现PLS
    ,可见这个算法的意义。
    if nargout <= 2
        [Xloadings,Yloadings] = simpls(X0,Y0,ncomp);
        
    elseif nargout <= 4
        [Xloadings,Yloadings,Xscores,Yscores] = simpls(X0,Y0,ncomp);
    ....
    end
    
    
    %------------------------------------------------------------------------------
    %SIMPLS Basic SIMPLS.  Performs no error checking.
    function [Xloadings,Yloadings,Xscores,Yscores,Weights] = simpls(X0,Y0,ncomp)
    
    [n,dx] = size(X0);
    dy = size(Y0,2);
    
    % An orthonormal basis for the span of the X loadings, to make the successive
    % deflation X0'*Y0 simple - each new basis vector can be removed from Cov
    % separately.
    
    %计算协方差矩阵
    
    Cov = X0'*Y0;
    for i = 1:ncomp
        % Find unit length ti=X0*ri and ui=Y0*ci whose covariance, ri'*X0'*Y0*ci, is
        % jointly maximized, subject to ti'*tj=0 for j=1:(i-1).
        %计算协方差的主成分方向   
        [ri,si,ci] = svd(Cov,'econ');
        ri = ri(:,1); ci = ci(:,1); si = si(1);
        %计算得分
        ti = X0*ri;
        normti = norm(ti); ti = ti ./ normti; % ti'*ti == 1
        %计算载荷
        Xloadings(:,i) = X0'*ti;
        
        qi = si*ci/normti; % = Y0'*ti
        Yloadings(:,i) = qi;
        
        if nargout > 2
            Xscores(:,i) = ti;
            Yscores(:,i) = Y0*qi; % = Y0*(Y0'*ti), and proportional to Y0*ci
            if nargout > 4
                Weights(:,i) = ri ./ normti; % rescaled to make ri'*X0'*X0*ri == ti'*ti == 1
            end
        end
    
        % Update the orthonormal basis with modified Gram Schmidt (more stable),
        % repeated twice (ditto).
        vi = Xloadings(:,i);
        for repeat = 1:2
            for j = 1:i-1
                vj = V(:,j);
                vi = vi - (vj'*vi)*vj;
            end
        end
        vi = vi ./ norm(vi);
        V(:,i) = vi;
    
        % Deflate Cov, i.e. project onto the ortho-complement of the X loadings.
        % First remove projections along the current basis vector, then remove any
        % component along previous basis vectors that's crept in as noise from
        % previous deflations.
    %更新协方差矩阵
        Cov = Cov - vi*(vi'*Cov);
        Vi = V(:,1:i);
        Cov = Cov - Vi*(Vi'*Cov);
    end
    
    if nargout > 2
        % By convention, orthogonalize the Y scores w.r.t. the preceding Xscores,
        % i.e. XSCORES'*YSCORES will be lower triangular.  This gives, in effect, only
        % the "new" contribution to the Y scores for each PLS component.  It is also
        % consistent with the PLS-1/PLS-2 algorithms, where the Y scores are computed
        % as linear combinations of a successively-deflated Y0.  Use modified
        % Gram-Schmidt, repeated twice.
        for i = 1:ncomp
            ui = Yscores(:,i);
            for repeat = 1:2
                for j = 1:i-1
                    tj = Xscores(:,j);
                    ui = ui - (tj'*ui)*tj;
                end
            end
            Yscores(:,i) = ui;
        end
    end
    
    
    %------------------------------------------------------------------------------
    %PLSCV Efficient cross-validation for X and Y mean squared error in PLS.
    function mse = plscv(X,Y,ncomp,cvp,mcreps,ParOptions)
    
    [n,dx] = size(X);
    
    % Return error for as many components as asked for; some columns may be NaN
    % if ncomp is too large for CV.
    mse = NaN(2,ncomp+1);
    
    % The CV training sets are smaller than the full data; may not be able to fit as
    % many PLS components.  Do the best we can.
    if isa(cvp,'cvpartition')
        cvpType = 'partition';
        maxncomp = min(min(cvp.TrainSize)-1,dx);
        nTest = sum(cvp.TestSize);
    else
        cvpType = 'Kfold';
    %    maxncomp = min(min( floor((n*(cvp-1)/cvp)-1), dx));
        maxncomp = min( floor((n*(cvp-1)/cvp)-1), dx);
        nTest = n;
    end
    if ncomp > maxncomp
        warning(message('stats:plsregress:MaxComponentsCV', maxncomp));
        ncomp = maxncomp;
    end
    
    % Cross-validate sum of squared errors for models with 1:ncomp components,
    % simultaneously.  Sum the SSEs over CV sets, and compute the mean squared
    % error
    CVfun = @(Xtr,Ytr,Xtst,Ytst) sseCV(Xtr,Ytr,Xtst,Ytst,ncomp);
    sumsqerr = crossval(CVfun,X,Y,cvpType,cvp,'mcreps',mcreps,'options',ParOptions);
    mse(:,1:ncomp+1) = reshape(sum(sumsqerr,1)/(nTest*mcreps), [2,ncomp+1]);
    
    
    %------------------------------------------------------------------------------
    %SSECV Sum of squared errors for cross-validation
    function sumsqerr = sseCV(Xtrain,Ytrain,Xtest,Ytest,ncomp)
    
    XmeanTrain = mean(Xtrain);
    YmeanTrain = mean(Ytrain);
    X0train = bsxfun(@minus, Xtrain, XmeanTrain);
    Y0train = bsxfun(@minus, Ytrain, YmeanTrain);
    
    % Get and center the test data
    X0test = bsxfun(@minus, Xtest, XmeanTrain);
    Y0test = bsxfun(@minus, Ytest, YmeanTrain);
    
    % Fit the full model, models with 1:(ncomp-1) components are nested within
    [Xloadings,Yloadings,~,~,Weights] = simpls(X0train,Y0train,ncomp);
    XscoresTest = X0test * Weights;
    
    % Return error for as many components as the asked for.
    outClass = superiorfloat(Xtrain,Ytrain);
    sumsqerr = zeros(2,ncomp+1,outClass); % this will get reshaped to a row by CROSSVAL
    
    % Sum of squared errors for the null model
    sumsqerr(1,1) = sum(sum(abs(X0test).^2, 2));
    sumsqerr(2,1) = sum(sum(abs(Y0test).^2, 2));
    
    % Compute sum of squared errors for models with 1:ncomp components
    for i = 1:ncomp
        X0reconstructed = XscoresTest(:,1:i) * Xloadings(:,1:i)';
        sumsqerr(1,i+1) = sum(sum(abs(X0test - X0reconstructed).^2, 2));
    
        Y0reconstructed = XscoresTest(:,1:i) * Yloadings(:,1:i)';
        sumsqerr(2,i+1) = sum(sum(abs(Y0test - Y0reconstructed).^2, 2));
    end
    
    

     

    展开全文
  • 最小二乘法MATLAB程序

    2009-10-11 11:05:27
    偏最小二乘法 最小二乘法是一种数学优化技术,它通过最小化误差的平方和找到一组数据的最佳函数匹配。 用最简的方法求得一些绝对不可知的真值,而令误差平方之和为最小。 通常用于曲线拟合。很多其他的优化问题也可...
  • 偏最小二乘法资料和matlab程序

    热门讨论 2009-04-03 21:47:54
    偏最小二乘法资料和matlab程序,介绍非常详细,对搞统计分析的有较大帮助!
  • 偏最小二乘回归方法用数学软件matlab进行开发,实现,这里是其源代码,还有一个例子。
  • 软件包,它同时使用偏最小二乘回归 (PLSR) 和自适应网络模糊推理系统 (ANFIS) 在一组自变量 ( X ) 和因变量之间建立预测模型。变量 ( Y )。 该模型采用如下两阶段方法: 阶段1: 第 1 部分:使用 PLSR,建立预测模型...
  • 这是一个介绍偏最小二乘法,希望对大家有所帮助
  • matlab偏最小二乘法及其检验

    万次阅读 多人点赞 2016-02-26 00:59:24
    偏最小二乘法,它通过最小化误差的平方和找到一组数据的最佳函数匹配。 用最简的方法求得一些绝对不可知的真值,而令误差平方之和为最小。 很多其他的优化问题也可通过最小化能量或最大化熵用最小二乘形式表达。...

    偏最小二乘法,它通过最小化误差的平方和找到一组数据的最佳函数匹配。 用最简的方法求得一些绝对不可知的真值,而令误差平方之和为最小。
    很多其他的优化问题也可通过最小化能量或最大化熵用最小二乘形式表达。

    • 偏最小二乘法
    clc % 清屏
    clear all; % 删除workplace变量
    close all; % 关掉显示图形窗口
    format short
    pz=[191 36  50  5   162 60
    189 37  52  2   110 60
    193 38  58  12  101 101
    162 35  62  12  105 37
    189 35  46  13  155 58
    182 36  56  4   101 42
    211 38  56  8   101 38
    167 34  60  6   125 40
    176 31  74  15  200 40
    154 33  56  17  251 250
    169 34  50  17  120 38
    166 33  52  13  210 115
    154 34  64  14  215 105
    247 46  50  1   50  50
    193 36  46  6   70  31
    202 37  62  12  210 120
    176 37  54  4   60  25
    157 32  52  11  230 80
    156 33  54  15  225 73
    138 33  68  2   110 43];%每一列为一个指标,每一行为一个样本
    
    mu=mean(pz); %求均值
    sig=std(pz); %求标准差
    rr=corrcoef(pz); %求相关系数矩阵
    data=zscore(pz); %数据标准化
    n=3; % n 是自变量的个数
    m=3; % m 是因变量的个数
    x0=pz(:,1:n);y0=pz(:,n+1:end);%定义自变量为前n列,因变量为n+1到m列
    e0=data(:,1:n);f0=data(:,n+1:end);%e0为自变量归一化值,f0为因变量归一化值
    num=size(e0,1);%求样本点的个数
    chg=eye(n); % w 到 w* 变换矩阵的初始化
    for i=1:n
        %计算 w,w* 和t 的得分向量,
        matrix=e0'*f0*f0'*e0;
        [vec,val]=eig(matrix); %求特征值和特征向量
        val=diag(val); %提出对角线元素,即特征值
        [val,ind]=sort(val,'descend');%降序排列,ind为排序后原下标序号
        w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量
        w_star(:,i)=chg*w(:,i); %计算w*的取值
        t(:,i)=e0*w(:,i); %计算成分ti 的得分
        alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算alpha_i
        chg=chg*(eye(n)-w(:,i)*alpha'); %计算w 到w*的变换矩阵
        e=e0-t(:,i)*alpha'; %计算残差矩阵
        e0=e;%将残差矩阵带进下次循环
        %计算ss(i)的值
        beta=[t(:,1:i),ones(num,1)]\f0; %求回归方程的系数
        beta(end,:)=[]; %删除回归分析的常数项
        cancha=f0-t(:,1:i)*beta; %求残差矩阵
        ss(i)=sum(sum(cancha.^2)); %求误差平方和
        %计算p(i)
        for j=1:num
            t1=t(:,1:i);f1=f0;
            she_t=t1(j,:);she_f=f1(j,:); %把舍去的第j 个样本点保存起来
            t1(j,:)=[];f1(j,:)=[]; %删除第j 个观测值
            beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数
            beta1(end,:)=[]; %删除回归分析的常数项
            cancha=she_f-she_t*beta1; %求残差向量
            p_i(j)=sum(cancha.^2);
        end
        p(i)=sum(p_i);
        if i>1
            Q_h2(i)=1-p(i)/ss(i-1);
        else
            Q_h2(1)=1;
        end
        if Q_h2(i)<0.0975
            fprintf('提出的成分个数r=%d',i);
            fprintf('   ');
            fprintf('交叉的有效性=%f',Q_h2(i));
            r=i;
            break
        end
    end
    beta_z=[t(:,1:r),ones(num,1)]\f0; %求Y 关于t 的回归系数
    beta_z(end,:)=[]; %删除常数项
    xishu=w_star(:,1:r)*beta_z; %求Y 关于X 的回归系数,且是针对标准数据的回归系数,每一列是一个回归方程
    mu_x=mu(1:n);mu_y=mu(n+1:end);
    sig_x=sig(1:n);sig_y=sig(n+1:end);
    for i=1:m
        ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数项
    end
    for i=1:m
        xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数,每一列是一个回归方程
    end
    sol=[ch0;xish] %显示回归方程的系数,每一列是一个方程,每一列的第一个数是常数项,每一列为一个因变量与自变量们的回归方程
    %此为还原为原始变量后的方程
    %% 感觉用途不大,用到的时候再查询怎么使用
    save mydata x0 y0 num xishu ch0 xish
    w1=w(:,1)
    w2=w(:,2)
    wx1=w_star(:,1)
    wx2=w_star(:,2)
    tx1=t(:,1)'
    tx2=t(:,2)'
    beta_z %回归系数
    xishu%系数矩阵,即未还原原始变量的系数,每一列为一个因变量与自变量的回归方程
    
    %% 用法:分别计算出第四列第五列第六列和前三列的线性回归关系,给出系数,系数以列的方式给出,
    %%sol分别为常数项系数,x1 x2 x3的系数
    
    • 生成mydata数据,进行偏最小二乘检验
    clc % 清屏
    clear all; % 删除workplace变量
    close all; % 关掉显示图形窗口
    format short
    load('mydata.mat')%mydata为计算偏最小二乘保存的数据集,可以用于检验
    %% 更直观的解释各个自变量的作用
    figure
    bar(xishu')%分别画出三个自变量对三个因变量标准化后回归方程的系数的的长度图
    axis tight
    hold on
    annotation('textbox',[0.26 0.14 0.086 0.07],'String',{'单杠'},'FitBoxToText','off');
    annotation('textbox',[0.56 0.14 0.086 0.07],'String',{'弯曲'},'FitBoxToText','off');
    annotation('textbox',[0.76 0.14 0.086 0.07],'String',{'跳高'},'FitBoxToText','off');%在指定位置加注释
    %% 拟合效果的确定
    %所有点都在对角线附近均匀分布,则效果较好
    ch0=repmat(ch0,num,1);%repmat起复制矩阵组合为新矩阵的作用
    yhat=ch0+x0*xish; %计算y 的预测值
    y1max=max(yhat);
    y2max=max(y0);
    ymax=max([y1max;y2max]);
    cancha=yhat-y0; %计算残差
    figure
    subplot(2,2,1)
    plot(0:ymax(1),0:ymax(1),yhat(:,1),y0(:,1),'*')
    title('单杠成绩预测')
    subplot(2,2,2)
    plot(0:ymax(2),0:ymax(2),yhat(:,2),y0(:,2),'O')
    title('弯曲成绩预测')
    subplot(2,1,2)
    plot(0:ymax(3),0:ymax(3),yhat(:,3),y0(:,3),'H')
    title('跳高成绩预测')
    %% 绘制拟合效果图和权重比重图
    
    • 结果:
      这里写图片描述
      这里写图片描述

    偏最小二乘回归≈多元线性回归分析+典型相关分析+主成分分析

    展开全文
  • matlab开发-偏最小二乘法和判别分析法。使用PLS进行判别分析的教程和工具。
  • pretreat.m,pretreat.m ,opls.m ,oscfearn.m loscwold.m ks.m pls.m Idapinv.m plslda.m lecr.m plscv.m, plsidacv.m lplscv.m, plsldacv.m, ecrcv.m plsdcv.m, plsldadcv.m plsmccv.m, plsldamccv.m,mcs.m ...
  • Matlab偏最小二乘法用于判,,,别分析
  • 偏最小二乘法代码偏最小二乘法 (PLS)、基于核的潜在结构正交投影 (K-OPLS) 和基于 NIPALS 的 OPLS 基于Yi Cao实现的PLS回归算法: 基于 的 K-OPLS 回归算法。 基于 R 包的 OPLS 实现,使用 NIPALS 分解循环。 安装 $...
  • 许多matlab最小二乘法的源程序,大家仔细看看m文件里的说明就可以了
  • matlab编写实现最小二乘法多元线性拟合,可以得到最终拟合方程,并画出预测的回归系数直方图
  • matlab偏最小二乘法

    千次阅读 2016-12-12 23:43:30
    偏最小二乘是建立表到表的线性拟合关系,然后预测的方法(处理高维数据),比如在光谱分析中,X是某物质的光谱样本构成的训练集,Y是对应的成分数据,x是要预测成分的光谱数据function [y5,S] = pls2(X,Y,x) ...

    偏最小二乘是建立表到表的线性拟合关系,然后预测的方法(处理高维数据),比如在光谱分析中,X是某物质的光谱样本构成的训练集,Y是对应的成分数据,x是要预测成分的光谱数据

    function [y5,S] = pls2(X,Y,x)
    %[y5,S] = pls2(X,Y,x)【参考自司守奎教材 偏最小而成回归章节代码】
    % X、Y 为模型生成样本,x为需要预测对象,y5是预测结果
    % S.sol为模型的矩阵
    S.name = char({'系数矩阵:S(第一行为常数项)';
    '相对误差:relativeError';    '均方根误差:RootMSE';
    '标准误差:StandardError'; '拟合优度:R2'
    });
    pz = [X,Y];%通过这里注意数据组织形式,一行为一个数据
    mu = mean(pz);sig = std(pz);
    mu=mean(pz);sig=std(pz); %求均值和标准差
    rr=corrcoef(pz); %求相关系数矩阵
    data=zscore(pz); %数据标准化
    n=size(X,2);m=size(Y,2); %n 是自变量的个数,m 是因变量的个数
    x0=pz(:,1:n);y0=pz(:,n+1:end);
    e0=data(:,1:n);f0=data(:,n+1:end);
    num=size(e0,1);%求样本点的个数
    chg=eye(n); %w 到 w*变换矩阵的初始化
    for i=1:n
        %以下计算 w,w*和 t 的得分向量,
        matrix=e0'*f0*f0'*e0;
        [vec,val]=eig(matrix); %求特征值和特征向量
        val=diag(val); %提出对角线元素
        [val,ind]=sort(val,'descend');
        w(:,i)=vec(:,ind(1)); %提出最大特征值对应的特征向量
        w_star(:,i)=chg*w(:,i); %计算 w*的取值
        t(:,i)=e0*w(:,i); %计算成分 ti 的得分
        alpha=e0'*t(:,i)/(t(:,i)'*t(:,i)); %计算 alpha_i
        chg=chg*(eye(n)-w(:,i)*alpha'); %计算 w 到 w*的变换矩阵
        e=e0-t(:,i)*alpha'; %计算残差矩阵
        e0=e;
        %以下计算 ss(i)的值
        beta=[t(:,1:i),ones(num,1)]\f0; %求回归方程的系数
        beta(end,:)=[]; %删除回归分析的常数项
        cancha=f0-t(:,1:i)*beta; %求残差矩阵
        ss(i)=sum(sum(cancha.^2)); %求误差平方和
        %以下计算 press(i)
        for j=1:num
            t1=t(:,1:i);f1=f0;
            she_t=t1(j,:);she_f=f1(j,:); %把舍去的第 j 个样本点保存起来
            t1(j,:)=[];f1(j,:)=[]; %删除第 j 个观测值
            beta1=[t1,ones(num-1,1)]\f1; %求回归分析的系数
            beta1(end,:)=[]; %删除回归分析的常数项
            cancha=she_f-she_t*beta1; %求残差向量
            press_i(j)=sum(cancha.^2);
        end
        press(i)=sum(press_i);
        if i>1
            Q_h2(i)=1-press(i)/ss(i-1);
        else
            Q_h2(1)=1;
        end
        if Q_h2(i)<0.0975
            fprintf('提出的成分个数 r=%d',i);
            r=i;
            break
        end
    end
    beta_z=[t(:,1:r),ones(num,1)]\f0; %求 Y 关于 t 的回归系数
    beta_z(end,:)=[]; %删除常数项
    xishu=w_star(:,1:r)*beta_z; %求Y关于X的回归系数,且是针对标准数据的回归系数,每一列是一个回归方程
    mu_x=mu(1:n);mu_y=mu(n+1:end);
    sig_x=sig(1:n);sig_y=sig(n+1:end);
    for i=1:m
        ch0(i)=mu_y(i)-mu_x./sig_x*sig_y(i)*xishu(:,i); %计算原始数据的回归方程的常数项
    end
    for i=1:m
        xish(:,i)=xishu(:,i)./sig_x'*sig_y(i); %计算原始数据的回归方程的系数,每一列是一个回归方程
    end
    sol=[ch0;xish]; %显示回归方程的系数,每一列是一个方程,每一列的第一个数是常数项
    %% 用模型进行预测
    y5= x*xish;
    for i=1:size(y5,2)
        y5(:,i)= y5(:,i)+ch0(i);
    end
    %% 求模型的误差以及拟合优度
    yy=X*xish;y = Y;
    for i=1:size(yy,2)
        yy(:,i)= yy(:,i)+ch0(i);
    end
        e2 = abs((yy-y)./y);%求模型的相对误差
        S.relativeError = e2;
        for i=1:size(yy,2)
              c1(i) = sqrt(    sum( (yy(:,i)-y(:,i)).^2  ) ./ size(yy,1) );%求模型的均方根误差
              c2(i) = sqrt(     sum((yy(:,i)-mean(y(:,i))).^2)./size(yy,1)      );%标准误差
              [~,~,~,~,temp] = regress(yy(:,i),[ones(size(y,1),1),y(:,i)]);
              c3(i) = temp(1);%拟合优度
        end
        S.RootMSE = c1;S.StandardError = c2;S.R2 =  c3;
    end
    
    展开全文
  • 阻尼最小二乘法matlab代码版本控制说明 分支:测试/验证 验证通过ij = 11到33 与MATLAB比较的结果-同意 在main.f90中添加了二进制文件写入; dat使用了太多空间。 计算应力分量 全局计算tau_ij ^ F 使用非集中式...
  • 脚本可以直接运行以重现偏最小二乘回归 (PLSR) 分析和每个基因分析的相关性。 此外,还包含一个 R 脚本以使用 和 制作直方图。 分析 调用偏最小二乘回归 (PLSR) 和每个基因相关性分析。 运行: script_call_PLSR_and...
  • 偏最小二乘法应用程序在matlab中实现
  • 偏最小二乘法代码PartialLeastSquares_forClassification 用于数据分类的偏最小二乘算法的 Matlab 实现。 这些代码是基于以下论文实现的: Alin, A. (2009) “当对象数量远大于变量数量时 PLS 算法的比较”,统计...

空空如也

空空如也

1 2 3 4 5
收藏数 84
精华内容 33
关键字:

偏最小二乘法matlab

matlab 订阅