unpls.pas

神经网络用于X荧光分析的源码。输入元素峰的强度

end;

//------------------------------------------------------------------------------
//transform the property values
procedure XScale(const Scale : Byte; //data scale - raw, main, auto
                 const Nrow : Word;  //count of calibration samples
                 const Ncol : Word;  //count of spectra channels
                 var X : TSMatr; //X values, to scale
                 var Xbar,Xstd : TDVect //mean and std.dev scale values
                 );
var
  I,J : Word;
begin

  //transform the properties to mean-centered values
  for J := 1 to Ncol do
  begin
     Xbar[J] := 0.0; //init: if not scale !!!
     if (Scale <> 0) and (NRow <> 0) then
     begin
        for I := 1 to Nrow do
           Xbar[J] := Xbar[J] + X[I,J]; //sum
        Xbar[J] := Xbar[J] / NRow; //mean

        for I := 1 to Nrow do
           X[I,J] := X[I,J] - Xbar[J]; //centering
     end; // if

     //transform the properties to autoscaled values
     Xstd[J] := 1.0; //init: if not scale !!!
     if (Scale = 2) and (NRow > 1) then
     begin
        Xstd[J] := 0.0;
        for I := 1 to NRow do
           Xstd[J] := Xstd[J] + X[I,J] * X[I,J];
        Xstd[J] := Sqrt(Xstd[J] / (nrow-1));
        if (Xstd[J] = 0.0) then Xstd[J] := 1.0;
        for I := 1 to Nrow do
           X[I,J] := X[I,J] / Xstd[J];
     end;//if
  end;//do

end;

//------------------------------------------------------------------------------
//get property with responce and property correlation values
procedure GetPropRespCorr;
var I,J,K : Word;
begin
{  if Nrow > 1 then
  begin
    //compute and print the property correlation matrix
    //Property Correlation Matrix
    for I := 1 to Ncol do
    begin
       for J := 1 to Ncol do
       begin
          XCor[J,I] := 0.0;
          for K := 1 to Nrow do
             XCor[J,I] := XCor[J,I] + Xscl[K,I] * Xscl[K,J];
          XCor[J,I] := XCor[J,I] / (Nrow - 1);
       end;// do
    end;// do


    //get correlation of each input property with response
    //Correlation of Properties with Response
    for J := 1 to NCol do
    begin
       YXCor[J] := 0.0;
       for I := 1 to NRow do
          YXCor[J] := YXCor[J] + Yscl[I] * Xscl[I,J];
       YXCor[J] := YXCor[J] / (Nrow - 1);
    end;// do
  end;}
end;

//------------------------------------------------------------------------------
//nonlinear iterative PLS analysis
//     "nipals" finds partial least squares regression coefficients
//     via the classical nonlinear iterative PLS algorithm devised
//     by Herman Wold
//
//     literature reference:
//
//     S. Rannar, F. Lindgren, P. Geladi and S. Wold, "A PLS Kernel
//     Algorithm for Data Sets with many Variables and fewer Objects.
//     Part 1: Theory and Algorithms", Journal of Chemometrics, 8,
//     111-125 (1994)  [see pages 113-114 for NIPALS implementation]

procedure Nipals(const Nrow : Word;  //count of calibration samples
                 const Ncol : Word;  //count of spectra channels
                 const Y : TSVect;   //Y
                 const X : TSMatr; //X
                 var Nrank : Word;  //rank of matrix
                 var b : TDMatr;//pls regression coefficients for prediction
                 verbose : Boolean
                 );
var
  I,J,K,M : Word;
  norm : Double;
  w : TDMatr; //array[1..MaxCol,1..MaxRank] of Double; //PLS Property Weights
  q : TDVect; //array[1..MaxRank] of Double; //PLS Response loading vector
  p : TDMatr;//array[1..MaxCol,1..MaxRank] of Double; // PLS Property Loadings
  u,t : TDVect;
  wpw : TDMatr; //array[1..MaxCol,1..MaxRank] of Double;
  pw  : TDMatr; //Rnk_Rnk;
  Yref : TSVect;  //Y: 1..Nrow
  Xref : TSMatr;  //X: 1..NRow, 1..NCol

begin
  Inc(Progress);
 // fmMain.ProgressBarRun.Position := Progress;
  //set the rank
  VerifyRank(Nrow,Ncol,NRank);

  SetLength(w,NCol+1,NRank+1);
  SetLength(q,NRank+1);
  SetLength(p,NCol+1,NRank+1);
  SetLength(u,NRow+1); //loadings
  SetLength(t,NRow+1); //scores
  SetLength(wpw,NCol+1,NRank+1);
  SetLength(pw,NRank+1,NRank+1);

  SetLength(Yref, NRow + 1); //y
  SetLength(Xref, NRow +1,NCol + 1);//x



  for I := 1 to nrow do Yref[I] := Y[I];

  for J := 1 to ncol do
    for I := 1 to nrow do
        Xref[I,J] := X[I,J];

  //perform the NIPALS-PLS1 algorithm for successive ranks - fortlaufende PC's (AR)
  for K := 1 to NRank do
  begin
     //init score vector, by PLS1 Y -vector is proxy for u
     for I := 1 to nrow do
       u[I] := Yref[I];
//Step 1; w = (X) * u /|(X) * u| - calculation of loading weight vector - w
     norm := 0.0; //norm
     for J := 1 to Ncol do
     begin
        w[J,K] := 0.0; //sum
        for I := 1 to Nrow do
          w[J,K] := w[J,K] + u[I]*Xref[I,J];
        norm := norm + w[J,K] * w[J,K];
     end;

     norm := Sqrt(norm);
     if (norm = 0.0) then
     begin
        Nrank := K - 1; //not sucsessfull
        Break; //
     end;

     for J := 1 to ncol do //normalization
        w[J,K] := w[J,K] / norm;

//Step 2: Calculation of score vector t = X * w
     norm := 0.0;
     for I := 1 to nrow do
     begin
        t[I] := 0.0;
        for J := 1 to ncol do
           t[I] := t[I] + Xref[I,J] * w[J,K];
        norm := norm + t[I] * t[I]; //normalization for step 3
        tt[I,K] := t[I];
     end;

//Step 3: Calculate of loading vector q = ty/|ty|
     q[K] := 0.0;
     for I := 1 to nrow do
        q[K] := q[K] + t[I] * Yref[I];
     q[K] := q[K] / norm;

//Step 4: Calculate of X-loadings p = Xt/tt
     for J := 1 to ncol do
     begin
        p[J,K] := 0.0;
        for I := 1 to nrow do
           p[J,K] := p[J,K] + t[I] * Xref[I,J];
        p[J,K] := p[J,K] / norm;
     end;

//Step 5: Deflation or updating X(k+1) = X(k) - tp; y(k+1) = y(k) - qt
     for J := 1 to ncol do
        for I := 1 to nrow do
           Xref[I,J] := Xref[I,J] - t[I] * p[J,K];

     for I := 1 to nrow do
        Yref[I] := Yref[I] - q[K] * t[I];

  end; //for K

  //Hier soll noch ein Schritt f黵 bestimmen von PC nummer !!!

//Step 6 or 7: Compute b0 and b for N PLS factors (PC's) to use then in te predictor
  // get the pls regression coefficients for increasing rank
  //pw
  for M := 1 to Nrank do //all PC's ???
  begin
    for I := 1 to M do
    begin
      for J := 1 to M do
      begin
         pw[J,I] := 0.0;
         for K := 1 to ncol do
            pw[J,I] := pw[J,I] + p[K,J] * w[K,I]
      end;
    end;

    Invert(M,pw);

    for I := 1 to ncol do
    begin
      for J := 1 to M do
      begin
        wpw[I,J] := 0.0;
        for K := 1 to M do
          wpw[I,J] := wpw[I,J] + w[I,K] * pw[K,J];
      end;
    end;

    for I := 1 to ncol do
    begin
      b[I,M] := 0.0;
      for J := 1 to M do
        b[I,M] := b[I,M] + wpw[I,J] * q[J];
    end; // for I
  end; //for M
end;

//------------------------------------------------------------------------------
//calculate PLS-R fit values and residual calibration variance
procedure CalcPLSFitVAlues(const Nrank : Word;
                           const Nrow : Word;  //count of calibration samples
                           const Ncol : Word  //count of spectra channels, column
                           );
var
  Yfit  : TSVect; //Y - fit value
  I,J,M : Word;
begin
  SetLength(Yfit,NRow + 1);
  //compute the PLSR model fit for each rank
  for M := 1 to NRank do
  begin
     for I := 1 to Nrow do
     begin
        Yfit[I] := 0.0;
        for J := 1 to Ncol do
           Yfit[I] := Yfit[I] + b[J,M] * (Xraw[I,J] - Xbar[J]) / Xstd[J];
        Yfit[I] := Yfit[I] * Ystd + Ybar; //recalculate to raw
        YY_fit[I,M] := Yfit[I]; //save to global fit var.
     end;

//     Rsq := Correlation(nrow,Yraw,Yfit,true); //correlarion coef. r
//     Rrms := RMSE(Nrow,Yraw,Yfit); //RMSEC
     ResCalVar[M] := ResVar(Nrow,Yraw,Yfit);//residual calibration variance
  end;
end;


//------------------------------------------------------------------------------
//get predictions from PLS model
procedure LOOPlSModel( const Scale : Byte; //data scale
                       const Nrow : Word;  //count of calibration samples
                       const Ncol : Word;  //count of spectra channels, column
                       var   Nrank : Word
                       );

var
  I,J,K,M,extrap : Word;
  nqrank,npred : Word;
  rsq,qsq,rave,qave,rrms,qrms,yqbar,yqstd : Double;
  Xqbar,Xqstd : TDVect; //scaling (mean,std) factors for LOOCV
  Yfit  : TSVect; //Y - fit value
  Ypred : TSVect; //Y - predicted value LOOCV
  //work for LOOCV
  Yq : TSVect;
  Xq : TSMatr;
  bq : TDMatr;
  Nameq : array of string;
begin
  Nqrank := NRank; //AR, init
//  SetLength(Yfit,NRow + 1);
  SetLength(Ypred,NRow + 1);
  SetLength(Xqbar,NCol + 1);//scaling (mean) for LOOCV
  SetLength(Xqstd,NCol + 1);//scaling (std) for LOOCV
  SetLength(Yq,NRow + 1);
  SetLength(Xq, NRow +1,NCol + 1);//x
  SetLength(bq, NCol + 1, MaxRank + 1);
  SetLength(Nameq,NRow + 1);

  //set the cross-validation to use the leave-one-out method
  Npred := Nrow - 1; //number of predicted in LOOCV

  for M := 1 to NRank do
  begin
  {
     for I := 1 to Nrow do
     begin
        Yfit[I] := 0.0;
        for J := 1 to Ncol do
           Yfit[I] := Yfit[I] + b[J,M] * X[I,J]; //Y-fit value
        Yfit[I] := Yfit[I] * Ystd + Ybar;        //recalculate to raw
        YY_fit[I,M] := Yfit[I]; //save to global fit var.
     end;

     Rsq := Correlation(nrow,Yraw,Yfit,true); //correlarion coef. r
     Rrms := RMSE(Nrow,Yraw,Yfit); //RMSEC
     ResCalVar[M] := ResVar(Nrow,Yraw,Yfit);//residual calibration variance
   }
     // get the cross-validated predictions for the current rank
     for I := 1 to Nrow do
     begin
        for J := 1 to I-1 do
        begin
           Nameq[J] := SmplLbl[J];
           Yq[J] := Yraw[J];
        end;

        for J := I to Npred do
        begin
           Nameq[J] := SmplLbl[J+1];
           Yq[J] := Yraw[J+1];  //J+1 shifted, LOO
        end;

        for K := 1 to ncol do
        begin