unitfp.pas

一个多元非线性回归分析源码以及其中的公式列表

      3:   // 差的方差最小 ,返回 差的方差
      begin
        dVP:=0;
        L:=high(dData)+1;
        setlength(dVvar,L);
        for i:=low(dData) to high(dData) do
        begin
          for j:=low(PDList) to high(PDList) do
            PDList[j].CAddress^:=dData[i,PDList[j].Index];

          Compile.GetValue(dV);
          dVP:=dVP+dV;
          dVvar[i]:=dV;
        end;

        dVP:=dVP/L;
        sum:=0;
        for i:=0 to L-1 do
        begin
          sum:=sum+sqr(dVvar[i]-dVP);
        end;
      end;
      4:   // 差的极差最小 ,返回 差的最大值减去最小值
      begin
        dVMax:=-MaxExtended;
        dVMin:=+MaxExtended;
        for i:=low(dData) to high(dData) do
        begin
          for j:=low(PDList) to high(PDList) do
            PDList[j].CAddress^:=dData[i,PDList[j].Index];

          Compile.GetValue(dV);
          if dV>dVMax then dVMax:=dV ;
          if dV<dVmin then dVMin:=dV;
        end;
        sum:=abs(dVMax-dVmin);
      end;

    end;
    //if (IsNan(Sum)) or (IsInfinite(Sum)) then
    //if (sum=Infinity) or (sum=NaN) or (sum=NegInfinity) then
    Temp:=sum;
    if ((PInt64(@Temp)^ and $7FF0000000000000) = $7FF0000000000000) then
      result:=sqrt(MaxExtended)
    else
      result:=sum;
  except
    result:=sqrt(MaxExtended);
  end;
end;

function getSubValue(const cList :array of TabcList;const pArray :array of extended;
                     const x:extended;const index0,index1:integer):extended;
var
  i,j   :integer;
  dV    :extended;
begin

  for i:=low(cList) to high(cList) do    //参数赋值
  begin
    cList[i].CAddress^:=pArray[i];
  end;


  for j:=low(PDList) to high(PDList) do
  begin
    if PDList[j].Index=index0 then PDList[j].CAddress^:=x;
    if PDList[j].Index=index1 then PDList[j].CAddress^:=0;
  end;

  Compile.GetValue(dV);
  result:=dV;

end;

procedure GetSubValue(const cList :array of TabcList;const pArray :array of extended;var dSubArray:array of extended);
var
  i,j   :integer;
  dV    :extended;
begin

  for i:=low(cList) to high(cList) do    //参数赋值
  begin
    cList[i].CAddress^:=pArray[i];
  end;
  for i:=low(dData) to high(dData) do
  begin
    for j:=low(PDList) to high(PDList) do
      PDList[j].CAddress^:=dData[i,PDList[j].Index];

    Compile.GetValue(dV);
    dSubArray[i]:=dV;
  end;
end;

function  strFind(const StrT1,StrT2:string;var index:integer):boolean; //在StrT1中是否有 标识符 StrT2
var
  i,L     :integer;
  Str0    :string;
  strX    :string;
begin
  Str0:=#0+uppercase(strT1)+#0+#0+#0;
  strX:=uppercase(strT2);
  L:=length(strX);
  for i:=2 to length(str0)-length(strX)-1 do
  begin
    if copy(str0,i,L)=strX then
    begin
      index:=i-1;
      result:=true;
      exit;
    end;
  end;
  index:=0;
  result:=false;
end;

function  GetKeyValue(const sKey:string):extended;
begin 
  result:=strtofloat(frmMain.ValueListEditorConfig.Values[sKey]);
end;


function GetDataAsStr(const sDataSet: TClientDataSet):string;
var
  strTemp : string;
  strT    : string;
  i       : integer;
  NullOk  : integer;

begin
  strTemp:='';
  sDataSet.First;
  while not (sDataSet.Eof) do
  begin
    strT:='';
    NullOk:=1;
    for i:=DataLength-1 downto 0 do
    begin
      if (NullOk=1) and (sDataSet.FieldByName('d'+inttostr(i)).Value=null)then
        continue
      else
      begin

        if NullOk=0 then
        begin
          if (sDataSet.FieldByName('d'+inttostr(i)).Value=null)then
            strT:=#9+strT
          else
            strT:=sDataSet.FieldByName('d'+inttostr(i)).AsString+#9+strT;
        end
        else
        begin
          NullOk:=0;
          if (sDataSet.FieldByName('d'+inttostr(i)).Value=null)then
            strT:=''
          else
            strT:=sDataSet.FieldByName('d'+inttostr(i)).AsString;
        end;
      end;
    end;
    strTemp:=strTemp+strT+#13+#10;
    sDataSet.Next;
  end;
  sDataSet.First;
  result:=strTemp;
end;

function GetStrFeildValue(var s:string):variant;
var
  i  :integer;

begin
  if (length(s)=0) or((length(s)>=2)and(s[1]=#13)and(s[2]=#10)) then
  begin
    result:=null;
    exit;
  end;

  if s[1]=#9 then  //tab key
  begin
    result:=null;
    s:=copy(s,2,length(s)-1);
  end
  else
  begin
    for i:=1 to length(s) do
    begin
      if (s[i]=#9) then
      begin
        result:=strtofloat(copy(s,1,i-1));
        s:=copy(s,i+1,length(s)-i);
        exit;
      end
      else if (s[i]=#13)then
      begin
        result:=strtofloat(copy(s,1,i-1));
        s:=copy(s,i,length(s)-i+1);
        exit;
      end;
    end;
    result:=strtofloat(s);
    s:='';
  end;    
end;


procedure SetDataAsStr(var sDataSet: TClientDataSet;const sData: string;const CellX:integer=0);
var
  str   :string;
  i     :integer;

  function GetCount(const s:string;const ins:char):integer;
  var
    j   :integer;
  begin
    result:=0;
    for j:=1 to length(s) do
    begin
      if s[j]=ins then inc(result);
    end;
  end;
begin
  str:=sData;
  sDataSet.First;

  while not (sDataSet.Eof) do
    sDataSet.Delete;
  for i:=1 to GetCount(str,#13) do
  begin
    sDataSet.Insert;
    sDataSet.Post;
  end;
  sDataSet.First;

  while true do
  begin
    if sDataSet.Eof then
      sDataSet.Insert
    else
      sDataSet.Edit;
    for i:=CellX to DataLength-1 do
        sDataSet.fieldByName('d'+inttostr(i)).Value:=GetStrFeildValue(str);
    sDataSet.Post;
    if length(str)=0 then exit;
    if (length(str)>=2)and(str[1]=#13)and(str[2]=#10) then
       str:=copy(str,3,length(str)-2);
    sDataSet.Next;
  end;


end;

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

procedure RunOptimize();
begin
  case ApWay of
    0:  RunOptimize0();  //随机爬山法
    1:  RunOptimize1();  //网格爬山法
    2:  RunOptimize2();  //最速下降法
    3:  RunOptimize3();  //最速下降网格爬山法
    4:  RunOptimize4();  //基因算法
  end;
end;

procedure RunOptimize0();   //随机爬山法
var
  MyTime      :extended;
  i,j         :integer;
  Grid        :int64;
  dTemp       :extended;
  MyabcList   :array of TabcList;
  MydMax      :extended;
  MyabcMax    :array of Extended;
  SumAgain    :extended;
  REAgainValue:extended;
  REAgain     :integer;
  properties0 :Extended;
  properties1 :Extended;
  properties2 :Extended;
begin
  //初始设定
  SetThreadPriority(GetCurrentThread(),THREAD_PRIORITY_LOWEST);
  randomize;
  MydMax:=MaxExtended;
  setlength(MyabcList,high(abcList)+1);
  setlength(abcMax,high(abcList)+1);
  setlength(MyabcMax,high(abcList)+1);
  for i:=low(MyabcList) to high(MyabcList) do
    MyabcList[i]:=abcList[i]; 
  ReAgain:=0;  
  REAgainValue:=-MaxExtended;
  properties0:=abs(GetKeyValue('随机爬山法.最小区域宽度'));
  properties1:=abs(GetKeyValue('随机爬山法.区域缩小倍数'));
  if properties1=0 then properties1:=1;
  properties2:=abs(GetKeyValue('随机爬山法.精度缩小倍数'));
  if properties2=0 then properties2:=1;
  //循环
  NowTime:=windows.GetTickCount();
  MyTime:=NowTime;
  while MyTime=NowTime do   //不按“停止”键就一直运行
  begin  
      application.ProcessMessages;
    Grid:=1;
    for i:=low(MyabcList) to high(MyabcList) do
      Grid:=Grid*MyabcList[i].Precision;
    for i:=0 to (Grid+1000) do
    begin
      for j:=low(MyabcList) to high(MyabcList) do
      begin
        if MyabcList[j].PType=2 then
          MyabcList[j].Default:=random(trunc(MyabcList[j].maxV-MyabcList[j].MinV)+1)+MyabcList[j].MinV
        else
          MyabcList[j].Default:=random*(MyabcList[j].maxV-MyabcList[j].MinV)+MyabcList[j].MinV;
      end;
      dTemp:=abs(GetFBValue(MyabcList));

      if dTemp<MydMax then    //局部最大
      begin
        MydMax:=dTEmp;
        for j:=low(MyabcMax) to high(MyabcMax) do
          MyabcMax[j]:=MyabcList[j].Default;
        application.ProcessMessages;
      end;
      if dTemp<dMax then     //全局最大
      begin
        dMax:=dTEmp;
        for j:=low(abcMax) to high(abcMax) do
          abcMax[j]:=MyabcList[j].Default;
        application.ProcessMessages;
      end;

      application.ProcessMessages;
      while runPause=true do
      begin
        application.ProcessMessages;
        sleep(1);
      end;
      application.ProcessMessages;
      if MyTime<>NowTime then exit;
      application.ProcessMessages;
    end;

    SumAgain:=0;
    for j:=low(MyabcList) to high(MyabcList) do
    begin
      if MyabcList[j].PType=2 then
      begin
        if abs(MyabcList[j].MaxV-MyabcList[j].MinV)>1 then
          sumAgain:=sumAgain+abs(MyabcList[j].MaxV-MyabcList[j].MinV);
      end
      else
        sumAgain:=sumAgain+abs(MyabcList[j].MaxV-MyabcList[j].MinV);
    end;
    if sumAgain=REAgainValue then
    begin
       REAgain:=REAgain+1;
    end
    else
    begin
      REAgainValue:=sumAgain;
      REAgain:=0;
    end;

    if (sumAgain<properties0) or (REAgain=10) then      //重来
    begin
      for i:=low(MyabcList) to high(MyabcList) do
        MyabcList[i]:=abcList[i];
      MydMax:=MaxExtended;
    end
    else
    begin
      for j:=low(MyabcList) to high(MyabcList) do
      begin
        MyabcList[j].MinV:=max(MyabcList[j].MinV,MyabcMax[j]-(MyabcList[j].MaxV-MyabcList[j].MinV)/(properties1)/2);//属性0
        MyabcList[j].MaxV:=min(MyabcList[j].maxV,MyabcMax[j]+(MyabcList[j].MaxV-MyabcList[j].MinV)/(properties1)/2);
        if MyabcList[j].PType=2 then
        begin
          MyabcList[j].MinV:=trunc(MyabcList[j].MinV);
          MyabcList[j].MaxV:=trunc(MyabcList[j].MaxV);
        end;
        MyabcList[j].Precision:=trunc(max(3.0,MyabcList[j].Precision/properties2));//属性1
      end;
    end;
    
  end;

end;

procedure RunOptimize1();   //网格爬山法
var
  MyTime      :extended;
  i,j         :integer;
  Grid        :int64;
  dTemp       :extended;
  MyabcList   :array of TabcList;
  MydMax      :extended;
  MyabcMax    :array of Extended;
  SumAgain    :extended;
  iAgain      :extended;
  ReAgain     :integer;