📄 unarmath.pas
字号:
unit unARMath;
//Writen by Dr. A.Ritter, 2002 - 2004
interface
uses Windows,Math;
type
TFloatMatrix = array of array of Double;
TFloatVector = array of Double;
TFitting_Formula = function(StartPosition, StopPosition : LongInt;
FittedParams, FittedValues : array of Double
) : Boolean;
function MaxFloatArray(
const StartPosition, StopPosition : LongInt; // > 0
var MaxPosition : LongInt;
const FloatArray : array of Double //start by 0
) : Double;
function MinFloatArray(const StartPosition, StopPosition : LongInt; // > 0
var MaxPosition : LongInt;
const FloatArray : array of Double //start by 0
) : Double;
function CalcWeigtArray(
const PoisonWeighted : Boolean;
const StartPosition, StopPosition : LongInt; // > 0
const MeasValues_Y : array of Double;
var WeightArray : array of Double
) : Boolean;
function Residuals( const PoisonWeighted : Boolean;
const StartPosition, StopPosition : LongInt; // > 0
const MeasValues_Y, FittedValues : array of Double;
var ResidsValues : array of Double
) : Boolean;
function Autocorrelation(
const PoisonWeighted : Boolean;
const StartPosition, StopPosition : LongInt; // > 0
const MeasValues_Y, FittedValues : array of Double;
var AutoCorrValues : array of Double;
var Count : LongInt
) : Boolean;
function Standard_Deviation(
const PoisonWeighted : Boolean;
const StartPosition, StopPosition : LongInt; // > 0
const FreedomDegree : LongInt;
const MeasValues_Y, FittedValues : array of Double;
const WeightArray : array of Double
) : Double;
function ResidsToDistribution(
const StartY,StopY : Double;
const StartPosition, StopPosition : LongInt; // > 0
const ChanCount : LongInt;
const ResidsValues : array of Double;
var Step : Double;
var DistrX, DistrY : array of Double
) : Boolean;
procedure Fitting(
const Fitting_Formula_No : LongInt;
const StartPosition, StopPosition : LongInt; // > 0
const PoisonWeighted : Boolean;
const IterStep : Double; //Iteration step
var StdDev : Double;
var Error_Count : Longint;
const Epsilon : Double; //Iteration limit
const MaxIter : LongInt; //Iteration limit
var Iter : LongInt;
const MeasValues_Y : array of Double;
const MeasValues_X : array of Double;
var FittedValues : array of Double;
const ParamCount : LongInt; //XX
var FittedParams : array of Double; //XX
var FittedParamsStdDev : array of Double
);
procedure CalcFitFunction(
const Fitting_Formula_No : LongInt;
const StartPosition, StopPosition : LongInt; // > 0
var FittedValues : array of Double;
const FittedParams : array of Double //XX
);
function LinRegression(const X,Y : array of double;
const PointCount : Word;
var Slope,Offset,Correlation,StdDev,Bias : Double)
: Boolean;
procedure AddParamsForCalibrFit(const RV,LV : TFloatMatrix;
const SC, RVC, LVC : Word);
//prepare NLMR
procedure NLMRprepare(const RV : TFloatMatrix; //row values
const PP : array of Boolean; //parameter is positiv
ParamCount, //Coefficients count
SamplesCount : Word);
implementation
var //espessially for calibration fitting, additionally params
RowValues,LabValues : array of array of Double;
SamplesCount,RowValueCount,LabValueCount : Word;
ParamIsPositiv : array of Boolean;
UsePosConstrains : Boolean;
//==============================================================================
// Fitting Functions
//==============================================================================
//------------------------------------------------------------------------------
//1 function Gauss
function Gauss(
const StartPosition, StopPosition : LongInt;
const FittedParams : array of Double;
var FittedValues : array of Double
) : Boolean;
var Coef : Double;
J : word;
W2,X2 : double;
begin
// FittedParams[0] = Centroid
// FittedParams[1] = Width
// FittedParams[2] = Area
Result := True;
try
for J := StartPosition to StopPosition do FittedValues[J] := 0.0;
if StopPosition <= StartPosition then Exit;
if FittedParams[1] <> 0 then
Coef := FittedParams[2] / (FittedParams[1] * Sqrt(PI / 2))
else
Exit; //0 lassen
W2 := FittedParams[1] * FittedParams[1];
for J := StartPosition to StopPosition do
begin
X2 := 2 * (J - FittedParams[0]) * (J - FittedParams[0]);
FittedValues[J] := Coef * Exp(- X2 / W2);
end;
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
//2 Polynom function, the power is dependet on the param. array dimension
function Polynom(
const StartPosition, StopPosition : LongInt;
const FittedParams : array of Double;
var FittedValues : array of Double
) : Boolean;
var P : Double;
ParamCount,J,K : LongInt;
begin
Result := True;
try
ParamCount := Length(FittedParams);
for J := StartPosition to StopPosition do
begin
P := 1.0;
FittedValues[J] := 0.0;
for K := 0 to ParamCount - 1 do
begin
FittedValues[J] := FittedValues[J] + FittedParams[K] * P;
P := P * J;
end;
end;
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
//22.08.2005
function Offset(
const StartPosition, StopPosition : LongInt;
const X : array of Double;
const FittedParams : array of Double; //ax + b; a,b
var FittedValues : array of Double
) : Boolean;
var J : LongInt;
begin
Result := True;
try
for J := StartPosition to StopPosition do
FittedValues[J] := FittedParams[0] + X[J];
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
function Linear(
const StartPosition, StopPosition : LongInt;
const X : array of Double;
const FittedParams : array of Double; //ax + b; a,b
var FittedValues : array of Double
) : Boolean;
var J : LongInt;
begin
Result := True;
try
for J := StartPosition to StopPosition do
FittedValues[J] := FittedParams[0] + FittedParams[1] * X[J];
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
//7
function Quadr(
const StartPosition, StopPosition : LongInt;
const X : array of Double;
const FittedParams : array of Double; //ax + b; a,b
var FittedValues : array of Double
): Boolean;
var J : LongInt;
begin
Result := True;
try
for J := StartPosition to StopPosition do
FittedValues[J] := FittedParams[0] * X[J] * X[J] + FittedParams[1] * X[J]
+ FittedParams[2];
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
{ de Jongh function for the concentrations calculation:
[D + E * R] * [1 + sum(Aj*Cj)]
ParamCount in Fitting call defined count of search/fitted parameters
for the case Alphas = fixed the ParamCount in Fitting call schall be 2
in another case this value is 2 + count of influence concentrations Cj
the fitted parameters are saved in FittedParams array as:
D at the 0 Position
E at the 1 Position
a1 at the 2 Position
...
}
function deJonghConc(
const FittedParams : array of Double;
var FittedValues : array of Double
) : Boolean;
var J,K : LongInt;
begin
Result := True;
try
for J := 0 to SamplesCount - 1 do //Samples
begin
FittedValues[J] := 0; //sum
for K := 0 to LabValueCount - 1 do
FittedValues[J] := FittedValues[J] + FittedParams[K+2] * LabValues[K,J]; // Cj , J := 1 to m concentrations
FittedValues[J] := (FittedParams[0] + FittedParams[1] * RowValues[0,J]) *
(1 + FittedValues[J]);
end;
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
function dJM2(
const FittedParams : array of Double;
var FittedValues : array of Double
) : Boolean;
var J,K : LongInt;
DC : Double;
AC : Double;
begin
Result := True;
try
for J := 0 to SamplesCount - 1 do //Samples
begin
DC := FittedParams[0] + FittedParams[1] * RowValues[0,J];
AC := FittedParams[2] * DC;
if AC <> 1 then AC := 1 - AC;
FittedValues[J] := 0; //sum
for K := 0 to LabValueCount - 1 do
FittedValues[J] := FittedValues[J] + FittedParams[K+3] * LabValues[K,J]; // Cj , J := 1 to m concentrations
FittedValues[J] := DC / AC * (1 + FittedValues[J]);
end;
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
//add params needed for calibration
procedure AddParamsForCalibrFit(const RV,LV : TFloatMatrix;
const SC, RVC, LVC : Word);
var I,J : Word;
begin
SamplesCount := SC;
RowValueCount := RVC;
LabValueCount := LVC;
SetLength(RowValues,RowValueCount,SamplesCount);
SetLength(LabValues,LabValueCount,SamplesCount);
for I := 0 to RowValueCount - 1 do
for J := 0 to SamplesCount - 1 do
RowValues[I,J] := RV[I,J];
for I := 0 to LabValueCount - 1 do
for J := 0 to SamplesCount - 1 do
LabValues[I,J] := LV[I,J];
end;
//------------------------------------------------------------------------------
//Peak Resolution in keV
function PeakResolFn(
const StartPosition, StopPosition : LongInt;
const X : array of Double;
const FittedParams : array of Double; //ax + b; a,b
var FittedValues : array of Double
) : Boolean;
const
NF : Double = 2.354820045; //2 * SQRT( 2 * ln(2))
EpsSi : Double = 3.85; //eV, Si Detector
var
J : LongInt;
begin
Result := True;
try
for J := StartPosition to StopPosition do
FittedValues[J] :=
SQRT(SQR(FittedParams[0] / NF) + (EpsSi / 1000) * FittedParams[1] * X[J]);
except
Result := False;
end;
end;
//------------------------------------------------------------------------------
// 8 non linear multiple regression fitted function
//------------------------------------------------------------------------------
//prepare NLMR
procedure NLMRprepare(const RV : TFloatMatrix; //row values
const PP : array of Boolean; //parameter is positiv
ParamCount, //Coefficients count
SamplesCount : Word);
var I,J : Word;
begin
SetLength(ParamIsPositiv,ParamCount);
SetLength(RowValues,ParamCount,SamplesCount);
for I := 0 to ParamCount - 1 do
begin
ParamIsPositiv[I] := PP[I];
if ParamIsPositiv[I] then UsePosConstrains := True;
end;
for I := 0 to ParamCount - 1 do
for J := 0 to SamplesCount - 1 do
RowValues[I,J] := RV[I,J];
end;
//------------------------------------------------------------------------------
//function itself
function NLMRcalc(const StartPosition, StopPosition : LongInt;
const ParamCount : LongInt;
var FittedParams: array of Double;
var FittedValues: array of Double)
: Boolean;
var I, J : Word;
begin
Result := True;
if UsePosConstrains then
begin
for I := 0 to ParamCount - 1 do
if ParamIsPositiv[I] then
FittedParams[I] := ABS(FittedParams[I]);
end;
try
//Records (Samples); StartPosition = 0; SamplesCount - 1 = StopPosition
for I := StartPosition to StopPosition do
begin
FittedValues[I] := 0; //sum
//Coefficients; CoeffCout = ParamCount
for J := 0 to ParamCount - 1 do
FittedValues[I] := FittedValues[I] + FittedParams[J] * RowValues[J,I];
end;
except
Result := False;
end;
end;
//==============================================================================
// Procedurs
//==============================================================================
//Distribution
function ResidsToDistribution;
var I,J : LongInt;
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -