unpls.pas

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

           for J := 1 to I-1 do
              Xq[J,K] := Xraw[J,K];

           for J := I to npred do
              Xq[J,K] := Xraw[J+1,K];
        end;

        Yscale (scale,npred,yq,yqbar,yqstd);
        Xscale (scale,npred,ncol,xq,xqbar,xqstd);
        Nipals (npred,ncol,yq,xq,nqrank,bq,false);

        //if nqrank <> nrank then ...
        // call simpls (npred,ncol,yq,xq,nqrank,bq,verbose)
        //berechnen Pred Wert und reskalieren falls notwendig (yqstd,xqstd <> 1,
        //x-,yqbar <> 0)
        Ypred[I] := 0.0;
        for J := 1 to Ncol do
          Ypred[I] := Ypred[I] + bq[J,M] *
            (Xraw[I,J] - Xqbar[J]) / Xqstd[J];
        Ypred[I] := Ypred[I] * Yqstd + Yqbar;
        yy_prd[I,M] := Ypred[I]; //save

        if Abort then
          if MessageDlg('Abort PLS-Regression?',mtConfirmation,[mbYes,mbNo],0) =
            mrYes then Exit;
     end; //for I

//     Qsq := Correlation(Nrow,Yraw,Ypred,True);
//     Qrms := RMSE(Nrow,Yraw,Ypred); //RMSEP !!!
     ResValVar[M] := ResVar(Nrow,Yraw,Ypred);

  end; //end for rank
end;

//------------------------------------------------------------------------------
//find the minimum and maximum value for each property
//to know if the predicted values are out or in the range
//last change 12.07.2005
procedure FindXRange(const NRow,NCol : Word;
                     const X : TSMatr);
var I,J : Word;
begin
  for I := 1 to Ncol do
  begin
    Xmax[I] := -1.0E6;
    Xmin[I] :=  1.0E6;
    for J := 1 to Nrow do
    begin
       Xmax[I] := Max(Xmax[I],X[J,I]);
       Xmin[I] := Min(Xmin[I],X[J,I]);
    end;
  end;
end;

//------------------------------------------------------------------------------
//procedure predicted PLS-R values Yp for Nrow samples
//04.06.2005
function PredictY( const Nrow : Word;  //count of predicted samples
                   const Ncol : Word;  //count of column in X variables matrix
                   const Nrnk : Word;  //prdivtion rank or PC
                   const Ybar, Ystd : Double; //data centering vor Y
                   const Xbar,Xstd : TDVect;  //data centering vor X
                   const Xraw  : TSMatr;  //measurements for predicted set, row
                   const b : TDMatr //coeff. for Prediction
                   ) : Double; //predicted values
var I,J : Word;
begin
  Result := 0.0;
  for J := 1 to Ncol do
    Result := Result + b[J,Nrnk] * //calculate
      (Xraw[Nrnk,J] - Xbar[J]) / Xstd[J]; //raw to center
  Result := Result * Ystd + Ybar;  //center to raw
end;



//------------------------------------------------------------------------------
//PLS Model validation with an external data set
//created Jun 2005
//last change 13.07.2005
procedure ExtPlsModel( const Rank : Word; //prediction rank
                       const NCol : Word; //data column number
                       const Nexv : Word //count of predicted values
                      );
var Extrap : Word; //extrapolation point counter
    Yprd : TSVect;
    M,K,J : Word;
    qsq, qrms : Double;
begin
  SetLength(Yprd,Nexv + 1);
  for M := 1 to Rank do // PCs
  begin
     // compute the PLSR model fit for the test set
     if (Nexv <> 0) then
     begin
        for K := 1 to Nexv do //samples
        begin
           Yprd[K] := 0.0;
           for J := 1 to Ncol do //calculate
             Yprd[K] := Yprd[K] + b[J,M] * (Xraw[M,J] - Xbar[J]) / Xstd[J]; //raw to center
           Yprd[K] := Yprd[K] * Ystd + Ybar; //rescale to row (init Ystd = 1, ybar = 0)
           YY_prd[K,M] := Yprd[K];

          { Extrap := 0;
           //kompare
           for J := 1 to Ncol do
           begin
              if (Xexv[K,J] > Xmax[J]) then Inc(extrap);
              if (Xexv[K,J] < Xmin[J]) then Inc(extrap);
           end;

           //falls X Werte nicht in Kalibration waren dann Extrapolation Warning
           if (Extrap = 0) then
           begin

           end
           else begin

           end;}
        end;
//        qsq := Correlation(Nexv,Yexv,Yfit,True);
//        qrms := RMSE(Nexv,Yexv,Yfit);
        ResValVar[M] := ResVar(Nexv,Yraw,Yprd);
     end;// if
  end;
end;

//------------------------------------------------------------------------------
//"invert" inverts a matrix using the Gauss-Jordan method
procedure Invert(const n : Word;
                 var   a : TDMatr);
var
  I,J,K : Word;
  icol,irow : Word;
  ipivot, indxc, indxr : array[1..MaxRank] of Integer;
  big,temp,pivot : Double;
begin

  //check to see if the matrix is too large to handle
  if (n > MaxRank) then
  begin
    MessageDlg('The matrix is too large to be inverted. Aborted',mtError,[mbOk],0);
    Exit;
  end;

  //perform matrix inversion via the Gauss-Jordan algorithm

  for I := 1 to n do
     ipivot[I] := 0;  //Drehpunkt - pivot

  for I := 1 to n do
  begin
     big := 0.0;
     for J := 1 to n do
     begin
        if (ipivot[J] <> 1) then
        begin
           for K := 1 to n do
           begin
              if (ipivot[K] = 0) then
              begin
                 if (abs(a[J,K]) >= big) then
                 begin
                    big := abs(a[J,K]);
                    irow := J;
                    icol := K;
                 end;
              end
              else begin
                 if (ipivot[K] > 1) then
                 begin
//                     write (*,20)
//   20                format (/,' INVERT  --  Cannot Invert',
//     &                          ' a Singular Matrix')
                    Exit; //????
                 end;
              end;
           end;
        end;// if
     end;
     ipivot[icol] := ipivot[icol] + 1;
     if (irow <> icol) then
     begin
        for J := 1 to n do
        begin
           temp := a[irow,J];
           a[irow,J] := a[icol,J];
           a[icol,J] := temp;
        end;
     end;// if
     indxr[I] := irow;
     indxc[I] := icol;
     if (a[icol,icol] = 0.0) then
     begin
//            write (*,30)
//   30       format (/,' INVERT  --  Cannot Invert a Singular Matrix')
        Exit;
     end;// if
     pivot := a[icol,icol];
     a[icol,icol] := 1.0;
     for J := 1 to n do
        a[icol,J] := a[icol,J] / pivot;

     for J := 1 to n do
     begin
        if (J <> icol) then
        begin
           temp := a[J,icol];
           a[J,icol] := 0.0;
           for K := 1 to n do
              a[J,K] := a[J,K] - a[icol,K] * temp;
        end;// if
     end;
  end;

  for I := n downto 1 do
  begin
     if (indxr[I] <> indxc[I]) then
     begin
        for K := 1 to n do
        begin
           temp := a[K,indxr[I]];
           a[K,indxr[I]] := a[K,indxc[I]];
           a[K,indxc[I]] := temp;
        end;
     end;// if
  end;
end;//invert


//------------------------------------------------------------------------------
//basic statistical measures
//------------------------------------------------------------------------------

function Mean(N : Word; X : TSVect) : Double;
var I : Word;
begin
  Result := 0.0;
  if N > 0 then
  begin
    for I := 1 to N do
      Result := Result + X[I];
    Result := Result / N
  end;
end;

//------------------------------------------------------------------------------

function Variance(N : Word; X : TSVect; FD : Boolean) : Double;
var I : Word;
    Xmean : Double;
begin
  Result := 0.0;
  if N > 0 then
  begin
    Xmean := Mean(N,X);
    for I := 1 to N do
      Result := Result + (X[I] - Xmean) * (X[I] - Xmean);
    if FD and (N > 1) then Result := Result / (N-1)
                      else Result := Result / N;
  end;
end;

//------------------------------------------------------------------------------

function StdDev(N : Word; X : TSVect; FD : Boolean) : Double;
begin
  Result := Sqrt(Variance(N, X, FD));
end;

//------------------------------------------------------------------------------

function Covariance(N : Word; X,Y : TSVect; FD : Boolean) : Double;
var I : Word;
    Ymean,Xmean : Double;
begin
  Result := 0;
  if N > 0 then
  begin
    Xmean := Mean(N,X);
    Ymean := Mean(N,Y);
    for I := 1 to N do
      Result := Result + (X[I] - Xmean) * (Y[I] - Ymean);
    if FD and (N > 1) then Result := Result / (N-1)
                      else Result := Result / N;
  end;
end;

//------------------------------------------------------------------------------
//AR, Correlation R (Pearson's correlation coefficient)
function Correlation(N : Word; X,Y : TSVect; FD : Boolean) : Double;
var Sx,Sy,CovXY : Double;
begin
  Result := 0;
  Sx := StdDev(N, X, FD);
  Sy := StdDev(N, Y, FD);
  CovXY := Covariance(N, X, Y, FD);
  if (Sx <> 0) and (Sy <> 0) then
    Result := CovXY / (Sx * Sy);
end;

//------------------------------------------------------------------------------
//Modeling Errors
//------------------------------------------------------------------------------
//Residual Variance
//valid for both calibration and prediction dependend on input values
function ResVar(N : Word; X,Y : TSVect) : Double;
var I : Word;
begin
  Result := 0.0;
  if N > 0 then
  begin
    for I := 1 to N do
      Result := Result + (X[I] - Y[I]) * (X[I] - Y[I]);
    Result := Result / N;
  end;
end;

//------------------------------------------------------------------------------
//AR RMSE error between given and calculated values
//valid for both RMSEC and RMSEP dependend on input values
function RMSE(N : Word; X,Y : TSVect) : Double;
begin
  Result := SQRT(ResVar(N,X,Y));
end;

//------------------------------------------------------------------------------

end.