uarrays.pas

采用Grubbs方式实现粗大误差的剔除.源码中包含一个自定义长度的数组.

unit Uarrays;

interface

uses
  SysUtils, Windows, Variants, Classes,math,UseLib,Dialogs;

type

  TDoubleArray = class(TPersistent)
  private
    FCapacity: Integer;
    FCount: Integer;
    FGrubbsCNT: Integer;
    FerrCNT: Integer;
    FOK: Boolean;
    function ReadOver:Boolean;
  public
    Data:Array of double;
    constructor Create(N: Integer);
    procedure Put(index: Integer; Value:Double);
    function  Get(index: Integer):Double;
   
    procedure Add(Value: Double);
    procedure Delete(Index: Integer);
    procedure Clear;
    
    property Count: Integer read FCount;
    property Capacity: Integer read FCapacity;
    property Over: Boolean read ReadOver;

    function AsString: String;
    function MaxIndex: Integer;
    function MinIndex: Integer;

    function DoGrbbs: double;
    Function DoGray:double;
    
    function gx: double;
    property GrubbsCNT: Integer read FGrubbsCNT;
    property errCNT: Integer read FerrCNT write FerrCNT;
    property OK: Boolean read FOK write FOK;
    function X: double;
    function Q: double;
  end;


implementation

const Grubbs5:array[3..30] of double
     =(1.15,1.46,1.67,1.82,1.94,2.03,2.11,2.18,2.23,2.28,2.33,2.37,2.41,2.44,
       2.48,2.50,2.53,2.56,2.58,2.60,2.62,2.64,2.66,2.74,2.81,2.87,2.91,2.96);

const Grubbs1:array[3..30] of double
     =(1.16,1.49,1.75,1.94,2.10,2.22,2.32,2.41,2.48,2.55,2.61,2.66,2.70,2.75,
       2.78,2.82,2.85,2.88,2.91,2.94,2.96,2.99,3.01,3.10,3.18,3.24,3.29,3.34);




constructor TDoubleArray.Create(N: Integer);
begin
  inherited Create;
  FCapacity := N;
  setLength(Data,N);
  FCount := 0;
  FGrubbsCNT:=0;
end;

procedure TDoubleArray.Put(index: Integer ; value: Double);
begin
 if Index < Capacity  then Data[Index]:=Value;
end;

function TDoubleArray.Get(index: Integer): Double;
begin
 result:=0;
 if Index < Capacity  then result:=Data[Index];
end;

procedure TDoubleArray.Add(Value: Double);
begin
if  Over then  delete(0);
     Data[Fcount]:=Value;
     inc(FCount);

{
 if  Over then begin
     Data[Fcount-1]:=Value;
     Fcount := Capacity ;
 end else begin
     Data[Fcount]:=Value;
     inc(FCount);
 end; 
 }
 
end;


procedure TDoubleArray.Clear;
 var i:integer;
begin
  FCount := 0; FGrubbsCNT:=0; FerrCNT:=0;  FOK:=False;
  for i:=0 to  Capacity -1 do Data[i]:=0;
end;

procedure TDoubleArray.Delete(Index: Integer);

begin
  if  (index< 0 ) or (Fcount=0)  or ( Index>=Fcount ) then exit;

  if  (index= High(Data) )then begin
      data[index]:= 0;
      dec(Fcount);
      exit;
  end;
  MoveMemory(@data[index], @data[index+1], (Capacity - index)* 8);
  
  data[Capacity-1]:=0;
  dec(Fcount);
 

end;

function TDoubleArray.AsString: String;
 var I: Integer;
begin
  Result :='';
  for I := Low(Data)  to High(Data) do   begin
    if I = High(Data)  then result:=Result + floattostr( Data[I] )
                       else result:=Result + floattostr( Data[I] )+#13#10;
  end;
end;


function TDoubleArray.ReadOver: Boolean;
begin
  result:=(Fcount>=Capacity);
end;


function TDoubleArray.MaxIndex: Integer;
var
  I: Integer; M:double;
begin
  result:=-1;
  if not Over then exit;

  Result := Low(Data);
  M := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do
    if M < Data[I] then begin
      M := Data[I];   Result := i;
    end;

end;

function TDoubleArray.MinIndex: Integer;
var
  I: Integer; M:double;
begin
  result:=-1;
  if not Over then exit;

  Result := Low(Data);
  M := Data[Low(Data)];

  for I := Low(Data) + 1 to High(Data) do
    if M > Data[I] then begin
      M := Data[I];   Result := i;
    end;
end;


function TDoubleArray.DoGrbbs: double;
 var x, q ,gx :double;
   N, MaxI,minI,I:Integer;
   old:string;
begin
 result:=-10000000000000;
 if not Over then exit;



 MaxI:=MaxIndex ;  MinI:=MinIndex ;
 x:=Mean(Data) ;   q:=StdDev(Data);
 if q=0 then exit;

 if (abs(Data[MaxI]-x)>abs(Data[MInI]-x)) then I:=MaxI else I:=minI;
 gx:= abs(Data[I]- x)/q;
 if gx>grubbs5[Capacity] then begin
  // showerr('Index:'+inttostr(I)+'  DATA[I]: '+ floattostr(Data[I]));
  result:=Data[I];
  Inc(FGrubbsCNT);
  delete(I);
 end;

end;


function TDoubleArray.gx: double;
 var x, q  :double;
   N, MaxI,minI,I:Integer;
   old:string;
begin
 result:=0;
 if not Over then exit;
 MaxI:=MaxIndex ;  MinI:=MinIndex ;
 x:=Mean(Data) ;   q:=StdDev(Data);
 if q=0 then exit;
 
 if (abs(Data[MaxI]-x)>abs(Data[MInI]-x)) then I:=MaxI else I:=minI;
 gx:= abs(Data[I]- x)/q;

end;


function TDoubleArray.Q: double;
begin
  result:=0;
  if Over then result:=StdDev(data)  ;
end;

function TDoubleArray.X: double;
begin
  result:=0;
  if Over then result:=Mean(data)  ;
end;






/////////////////////////////////////

  const k=3.75;
procedure SortTo(SrcA: array of double; var DstA: array of double);
var
  I, J: Integer;      T  :Double;
begin
  for I := High(SrcA) downto Low(SrcA) do DstA[i]:=SrcA[i];
  for I := High(DstA) downto Low(DstA) do begin

    for J := Low(DstA) to High(DstA) - 1 do
      if DstA[J] > DstA[J + 1] then
      begin
        T := DstA[J];
        DstA[J] := DstA[J + 1];
        DstA[J + 1] := T;
      end;

 end;
end;

Function GetP(N:integer):Integer;
begin
   if odd(N) then result:=(N+1) div 2
             else result:= N div 2;
end;

Procedure MakeSA(A:array of double; var SA:  array of double ) ;
 var  i:integer;  T:double;
begin

 T:=0;
 for i:=low(A)  to  high(A) do begin
      T:= T+A[i];
      SA[i]:=T;
 end;
end;


Function MaxIsErr(RA:array of double; X,maxp:double;p:integer):Boolean;
   var T:double;  N:integer;
begin
   N:=high(RA)+1;
   T:=X+k*maxp/(N-p);
   Result:=   (RA[N-1]< X) or( RA[N-1] > T)  ;
end;

Function MinIsErr(RA:array of double; X,maxp:double;p:integer):Boolean;
   var T:double;
begin

   T:=X-k*maxp/p;
   Result:=   (RA[0]> X) or( RA[0] < T)  ;
end;


Function TDoubleArray.DoGray:double;
 var sortA,SumA:array of double;
  N, p, i:integer; X, MaxP:double;
begin
 result:=-10000000000000;
 if not Over then exit;
 setlength(SortA,Count);  setlength(sumA,Count); 
 X:=Mean(Data) ;  p:=Getp(Count);
 SORTto(Data,sortA);
 makeSA(sortA,SumA);
 Maxp:=abs(SumA[p-1]- X* p);

 if MinIsErr(sortA,X,maxp,P) then begin
    //showmessage('最小值是异常数据') ;
    i:= MinIndex ;
    result:=data[i];
    delete(i);
  end;
  if MaxIsErr(sortA,X,maxp,P) then begin
   // showmessage('最大值是异常数据') ;
    i:= MaxIndex ;
    result:=data[i];
    delete(i);
  end;


end;





end.