偏最小而乘法程序能运行，结果不对！请高手帮忙看看呀

kang14

A＝[0.301 0.422 0.516 0.542 0.501 0.415 0.336 0.282 0.251 0.237 0.234 0.226 0.211 0.187 0.153 0.116 0.080 0.051 0.029 0.015 0.007;0.225 0.323 0.414 0.471 0.498 0.496 0.481 0.466 0.432 0.405 0.366 0.325 0.271 0.219 0.170 0.124 0.084 0.052 0.030 0.015 0.007;0.205 0.292 0.358 0.388 0.366 0.340 0.307 0.300 0.314 0.345 0.375 0.384 0.377 0.349 0.289 0.226 0.155 0.100 0.058 0.030 0.014;0.140 0.216 0.299 0.381 0.453 0.516 0.558 0.572 0.551 0.511 0.451 0.384 0.313 0.247 0.185 0.129 0.087 0.053 0.030 0.015 0.007;0.129 0.194 0.260 0.319 0.367 0.410 0.445 0.479 0.501 0.509 0.509 0.488 0.446 0.380 0.380 0.233 0.161 0.101 0.058 0.030 0.015;0.109 0.158 0.201 0.227 0.238 0.251 0.266 0.297 0.343 0.389 0.438 0.464 0.467 0.429 0.361 0.281 0.197 0.124 0.072 0.038 0.018 ]
C=
[0.6 0.2 0.2;0.4 0.5 0.2;0.4 0.2 0.4;0.2 0.7 0.2;0.2 0.5 0.4;0.2 0.2 0.5 ]
Aun=A
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; %[0.5 0.097 0.101;0.397 0.402 0.099;0.301 0.100 0.299;0.101 0.598 0.100;0.098 0.397 0.302;0.099 0.102 0.399];%是样品的浓度矩阵 随样品而变化
     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



X等于C才对

[ 本帖最后由 zhangnan3509 于 2007-7-4 14:15 编辑 ]