偏最小二乘法的程序结果不对

lzq8281

clear
A=input('请输入校正吸光度矩阵A=');
C=input('请输入校正浓度矩阵 C=');
for j=1:21                         %j为波长点数
      Aj=sum(A(:,j))/6;           %9由样品数决定
      Sj=std(A(:,j));
      for i=1:6                     %i为样品数
       A(i,j)=(A(i,j)-Aj)/Sj;   %中心化和标准化
      end
end
A
for k=1:3                              %k为组分数
    Ck=sum(C(:,k))/6;;                 %9由样品数决定
    Sk=std(C(:,k));
        for i=1:6                   %i为样品数
           C(i,k)=(C(i,k)-Ck)/Sk ;        %标准化和中心化
       end
end
C
    h=0;                            %输出处理后的吸光度矩阵A和浓度矩阵C
for i=1:1000
    h=h+1;
    r=C(:,2);
    w=(r'*A)/(r'*r);               %w和它的转置交换了
    w=w/norm(w);
    t=A*w'/(w*w');
    q=t'*C/(t'*t);                 %q和它的转置交换了
    q=q/norm(q);
    r=C*q'/(q*q');
    for i=1:1000
        t0=t;
        w=(r'*A)/(r'*r);
        w=w/norm(w);
        t=A*w'/(w*w');
        q=t'*C/(t'*t);
        q=q/norm(q);
        r=C*q'/(q*q');    %r没有和它的转置交换
        if norm(t-t0)<=0.00000001*norm(t0)
                break
        end
     end
     v=(t'*A)/(t'*t);
     vn=v';
     vn=v/norm(v);   
     t=t*norm(v);
     w=w*norm(v);
     b=(r'*t)/(t'*t);
     a=A-t*v;
     c=C-b*t*q;
     if  norm(a-A)<=0.0001%*norm(a) %norm(a-A)<=10.^-4 %
            break ;
     end
     if   norm(c-C)<=0.0001%*norm(c)%norm(c-C)<=10.^-4 %
            break;
     end
     A=a;
     C=c;
end
     v;
     q;
     k=h;
   Aun=input('请输入测试液的吸收矩阵Aun=');
   for j=1:21                                  %j为波长点数
        Sj=std(Aun(:,j));
        Aunj=sum(Aun(:,j))/6;                   %6由样品数决定
          for i=1:6                               %i为样品数
              Aunj1=(Aun(i,j)-Aunj);                    %中心化
              Aun(i,j)=(Aunj1)/Sj;                    %标准化
          end
    end
     Cun=t*b*q; %是样品的浓度矩阵 随样品而变化
     m=0;
   for i=1:1000
       m=m+1;
       t=Aun*w';
       Aun=Aun-t*v;
       Cun=Cun+b*t*q;   
       if m>k
          break
       end
   end
  X=Cun
  for i=1:3
      Xi=sum(X(:,i))/6
      Sx=std(X(:,i))
            for j=1:6
          X(j,i)=X(j,i)*Sx+Xi;
      end
  end

  X

[ 本帖最后由 eight 于 2007-8-20 20:18 编辑 ]

lzq8281

刚才发急了，有些问题没说明清楚，主要是这个程序不能得到我预期想要的结果

chino

现在解决了吗？