unarmath.pas

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

          MessageBeep(0);
          Exit;
       end;
  2 :  if not Polynom(StartPosition, StopPosition, FittedParams, FittedValues) then
       begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;

  3:   if not Linear(StartPosition, StopPosition, MeasValues_X, FittedParams, FittedValues) then
       begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;

  4:  if not PeakResolFn(StartPosition, StopPosition, MeasValues_X, FittedParams, FittedValues) then
      begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;
  5:  if not deJonghConc(FittedParams, FittedValues) then
      begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
      end;
   6: if not dJM2(FittedParams, FittedValues) then
      begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
      end;
   7 : if not Quadr(StartPosition, StopPosition, MeasValues_X, FittedParams, FittedValues) then
       begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;
   8 : if not NLMRcalc(StartPosition, StopPosition, ParamCount,FittedParams,FittedValues) then
       begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;
  9 :  if not Offset(StartPosition, StopPosition, MeasValues_X, FittedParams, FittedValues) then
       begin
          Inc(error_Count);
          MessageBeep(0);
          Exit;
       end;
  end;

  for I:= StartPosition to StopPosition  do
    WorkVect[I] := FittedValues[I] - MeasValues_Y[I];

  Result := Standard_Deviation( PoisonWeighted,StartPosition, StopPosition,
                                ParamCount,MeasValues_Y, FittedValues,WeightArray);
{
  if UsePosConstrains then
  begin
    for I := 0 to ParamCount - 1 do
      if ParamIsPositiv[I] then
        Result := Result - (FittedParams[I] /ABS(FittedParams[I]) - 1);
    //Falls Positiv passiert nicht, falls negativ erh鰄t sich um 2.
  end;
}
end;

//------------------------------------------------------------------------------
// Gradient, derivatives are numeric calculated
function Derivatives : LongInt;
var  TmpV : array of Double;
     DX   : Double; //Step
     J,I  : Word;
     Save : Double;
begin
   Result := 0;
   SetLength(TmpV,StopPosition + 1);
   try
     for J := 0 to ParamCount - 1 do
     begin
       //Step define
       if FittedParams[J] <> 0 then Dx := FittedParams[J] * 1E-5 else Dx := 1E-5;
       //save
       Save := FittedParams[J];
       for I := StartPosition to StopPosition do TmpV[I] := FittedValues[I];
       //set new shifted value and calculate
       FittedParams[J] := FittedParams[J] + DX;
       Fit_Function;
       //back
       FittedParams[J] := Save;
       //calculate part. derivatives
       for I := StartPosition to StopPosition do
       begin
         WorkArr[J,I] := (FittedValues[I] - TmpV[I]) / Dx;
         FittedValues[I] := TmpV[I];
       end;
     end;
   except
     Inc(Result);
   end;
end;

//------------------------------------------------------------------------------
// New Matrix
function prNew_Matrix(var A : TFloatMatrix;
                      var Diag, G : TFloatVector) : LongInt;
var I,J,K   : word;
    Sum     : extended;
begin
  Result := 0;
  if ParamCount = 0 then Exit;
  try
    for J := 0 to ParamCount - 1 do
    begin
       for K:= 0 to ParamCount - 1 do
       begin
          Sum := 0;
          for I:= StartPosition to StopPosition do
             Sum := Sum +  WeightArray[I] * WorkArr[J,I] * WorkArr[K,I];
          A[J,K] := Sum;
       end; {for K}
    end; {for J}

    for J := 0 to ParamCount - 1 do
    begin
       Diag[J] := A[J,J];
       if Diag[J] = 0.0 then Diag[J] := 1.0;
    end;

    for J:= 0 to ParamCount - 1 do
    begin
      Sum := 0;
      for I:= StartPosition to StopPosition  do
        Sum := Sum -  WeightArray[I] * WorkVect[I] * WorkArr[J,I];
      G[J] := Sum / Sqrt(Diag[J]);
    end;

    for J := 0 to ParamCount - 1 do
      for K := 0 to ParamCount -1 do
        A[J,K] := A[J,K] / ( Sqrt(Diag[J]) * Sqrt(Diag[K]) );

  except
    Inc(Result);
  end;
end; {procedure}

//------------------------------------------------------------------------------
// Solve System
procedure prSolveMatrix(const A : TFloatMatrix;
                        const G,Diag : TFloatVector;
                        var DeltaX : TFloatVector;
                        var C : TFloatMatrix; Step : Double);
var J,K : Word;
    B   : TFloatVector;
begin
   SetLength(B,ParamCount);
   for J := 0 to ParamCount - 1 do
     for K := 0 to ParamCount - 1 do
       if K <> J then C[J,K] := A[J,K]
                 else C[J,K] := A[J,K] * (Step + 1.0); // K = J

   for J := 0 to ParamCount - 1 do B[J]:= G[J];

   Gauss_Jordan(ParamCount,C,B);

   for J := 0 to ParamCount - 1 do
   begin
    if Diag[J] > 0 then
       DeltaX[J] := B[J] / SQRT(Diag[J])
    else begin
       DeltaX[J] := 0;
       Inc(Error_Count)
     end;
  end;

end; {procedure}

//------------------------------------------------------------------------------
//Calculate StdDev for parameters
procedure CalcFittedParamsStdDev;
var J : Word;
    Save : Double;
begin
  for J := 0 to ParamCount - 1 do
    FittedParamsStdDev[J] := Sqrt(Abs(C[J,J]));
  if StopPosition <> StartPosition then
  begin
    if not PoisonWeighted then
    begin
       if (StopPosition - StartPosition) > (ParamCount + 1) then
         Save := SQRT(StdDev / (ABS(StopPosition - StartPosition) - ParamCount -1))
       else
         Save := SQRT(StdDev /(ABS(StopPosition - StartPosition)));
       for J := 0 to ParamCount - 1 do
         FittedParamsStdDev[J] := Save * FittedParamsStdDev[J];
    end;
  end;  

end;

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

function fnEnd_Iter : Boolean;
var J : word;
   S,Sum  : Byte;
begin
  Result := False;

  if (Step > 1000) or (Step < 1E-10) then
  begin
    Result := true;
    Exit;
  end;

  Sum := 0;
  for J := 0 to ParamCount - 1 do
  begin
    if (abs(DeltaX[J]) / (Tau + abs(Save_XX[J])) < Epsilon)
       then s := 0
    else s := 1;
    Sum := Sum + s;
  end;

  if Sum = 0 then
  begin
    Result := true;
    Exit;
  end;

  if Save_StdDev <> 0 then
    if Abs(F - Save_StdDev) / Save_StdDev < Epsilon then
    begin
      Result := true;
      Exit;
    end;

  Save_StdDev := F;

  if (Iter >= MaxIter) then Result := true;

end;
//------------------------------------------------------------------------------
procedure SetDimensions;
begin
  SetLength(WorkArr,ParamCount,StopPosition + 1);
  SetLength(WeightArray,StopPosition + 1);
  SetLength(WorkVect,StopPosition + 1);
  SetLength(A,ParamCount,ParamCount);
  SetLength(C,ParamCount,ParamCount);
  SetLength(DeltaX,ParamCount);
  SetLength(Diag,ParamCount);
  SetLength(G,ParamCount);
  SetLength(Save_XX,ParamCount);
end;
//------------------------------------------------------------------------------

BEGIN {MAIN PROCEDURE}
  //Init
  Error_Count := 0;
  Iter := 0;
  Step := IterStep;
  //Dim
  SetDimensions;

  if not CalcWeigtArray(PoisonWeighted,StartPosition, StopPosition, MeasValues_Y,WeightArray)
    then Inc(Error_Count);

  //save
  Save_Step := Step;
  for J := 0 to ParamCount - 1 do Save_XX[J] := FittedParams[J];
  F := Fit_Function;
  Save_StdDev := F;

  repeat
    Inc(Iter);    //Iterations counter
    Derivatives;
    prNew_Matrix(A,Diag,G); //form matrix

    {Step I}
    Step := Save_Step / StepDevider;                          //reduce step
    prSolveMatrix(A,G,Diag,Deltax,C,Step);                         //calculate
    for J := 0 to ParamCount - 1 do FittedParams[J] := Save_XX[J] + DeltaX[J];//New values

    F2 := Fit_Function; //new StdDev
    if F2 <= F then  //if success than save and continue
    begin
      Save_Step := Step;
      for J := 0 to ParamCount - 1 do Save_XX[J] := FittedParams[J];
      F := F2;
      Continue;
    end; {if}

    {Step ii}
    Step := Save_Step / StepDevider;  //if not success reduce step
    prSolveMatrix(A,G,Diag,Deltax,C,Step);                         //calculate
    for J := 0 to ParamCount -1 do FittedParams[J] := Save_XX[J] + DeltaX[J];
    F1 := Fit_Function; //recalculate
    if (F1 <= F) and (F2 > F) then //compare
    begin
      Save_Step := Step; //if success than save and continue
      for J := 0 to ParamCount - 1 do Save_XX[J] := FittedParams[J];
      F := F1;
      Continue;
    end; {if}

    {Step iii}
    if (F1 > F) and (F2 > F) then //if both previous not success
    begin
      W := 0.0001;
      repeat
        Step := Save_Step * Exp(W * Ln(StepDevider)); //inc slowly step
        prSolveMatrix(A,G,Diag,Deltax,C,Step);        //calculate
        for J := 0 to ParamCount - 1 do FittedParams[J] := Save_XX[J] + DeltaX[J];
        F3 := Fit_Function;
        if F3 <= F then
        begin  //if success than save and continue
          Save_Step := Step;
          for J := 0 to ParamCount - 1 do Save_XX[J] := FittedParams[J]; //new values
          F := F3;
        end
          else W := W * 2; //
      until (F3 <= F) or (Step > 1000);
    end; {if}
  until fnEnd_Iter;

 //Calculate normalized standard deviation
  StdDev := Standard_Deviation( PoisonWeighted, StartPosition, StopPosition,
                                ParamCount, MeasValues_Y, FittedValues,
                                WeightArray);
  CalcFittedParamsStdDev;
end; {PROCEDURE MARQUARDT}

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

procedure CalcFitFunction;
begin
  case Fitting_Formula_No of
  1 :  Gauss(StartPosition, StopPosition, FittedParams, FittedValues);
  2 :  Polynom(StartPosition, StopPosition, FittedParams, FittedValues);
  end;
end;

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

function LinRegression(const X,Y : array of double;
                       const PointCount : Word;
                       var Slope,Offset,Correlation,StdDev,Bias : Double) : Boolean;
var AverX, AverY  : Extended;
    XXSum,YYSum,XYSum : Extended;
    K : word;
    Sigma2 : Double;
begin
  Result := False;

  if PointCount < 3 then  Exit;

  AverX := 0;
  for K := 0 to PointCount - 1 do AverX := AverX + X[K];
  AverX := AverX / PointCount;

  AverY := 0;
  for K := 0 to PointCount - 1 do AverY := AverY + Y[K];
  AverY := AverY / PointCount;

  XXSum := 0;
  YYSum := 0;
  XYSum := 0;

  for K := 0 to PointCount - 1 do
  begin
     XYSum := XYSum + (X[K] - AverX) * (Y[K] - AverY);
     XXSum := XXSum + (X[K] - AverX) * (X[K] - AverX);
     YYSum := YYSum + (Y[K] - AverY) * (Y[K] - AverY);
  end;

  if (XXSum = 0) or (YYSum = 0) then Exit;

// Koeffizienten

  Slope := XYSum / XXSum;
  Offset := AverY - Slope * AverX;

//  Korrelation
  Correlation := XYSum / SQRT(XXSum * YYSum);

//  Streuung

  Sigma2 := YYSum / (PointCount - 1);
  Sigma2 := (PointCount - 1) / (PointCount - 2) * Sigma2 *
        (1 - Correlation * Correlation);

  //CHI2
  StdDev := 0;
  for K := 0 to PointCount - 1 do
    StdDev := StdDev + SQR(Y[K] - X[K]);
  StdDev := SQRT(StdDev/(PointCount-1));

  //Bias
  Bias := 0;
  for K := 0 to PointCount - 1 do
    Bias := Bias + (X[K] - Y[K]);
  Bias := SQRT(ABS(Bias) / PointCount );

  Result := True;
end;


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

end.