unitlinearregressfunction.pas

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

unit UnitLinearRegressFunction;

interface
  uses Valedit, SysUtils, ExtCtrls, Types, Graphics, windows, Forms;

  type TXX=record
      n      :integer;         {数据组数n        }
      A      :double;          {线性回归方程系数A}
      B      :double;          {线性回归方程系数B}
      r      :double;          {相关系数r        }
      Sr     :double;          {标准误差Sr       }
      Sa     :double;          {A的误差Sa        }
      Sb     :double;          {B的误差Sb        }

      Px     :double;          {x的平均值        }
      Py     :double;          {y的平均值        }
      Hx     :double;          {x的和            }
      Hy     :double;          {y的和            }
      Hxy    :double;          {x*y的和          }
      Hxx    :double;          {x*x的和          }
      Hyy    :double;          {y*y的和          }

      xMin   :double;          {x的最小值xMin    }
      yMin   :double;          {x的最大值xMax    }
      xMax   :double;          {y的最小值yMin    }
      yMax   :double;          {y的最大值yMax    }
  end;
  var XX     :Txx;
  
  const Infinity:double=1.7e308; {无穷大}

  procedure getVLEData(const VLE :TValueListEditor;var Datax :array of double;var Datay :array of double);
  function  getDataMin(const Data:array of double;const n:integer):double;
  function  getDataMax(const Data:array of double;const n:integer):double;
  function  getDataSum(const Data:array of double;const n:integer):double;
  function  getDataProductSum(const Datax:array of double;const Datay:array of double;const n:integer):double;
  function  getDataSr(const Datax:array of double;const Datay:array of double;const n:integer;const A:double;const B:double):double;
  procedure Calculate(const Datax:array of double;const Datay:array of double;var xx:Txx);

  procedure DarwPic(var Img :TImage;const Datax:array of double;const Datay:array of double;const xx:Txx);
  procedure SetVLEData(var VLE :TValueListEditor;const xx :Txx);overload;
  procedure SetVLEData(var VLE :TValueListEditor);overload;
  procedure VLEInsertRow(var VLE:TValueListEditor;const str:string;const x:double);

  function  GetFileNamePath(const FileName:string):string;

implementation

procedure getVLEData(const VLE :TValueListEditor;var Datax :array of double;var Datay :array of double);
var i   :integer;
    n   :integer;
    s   :string;
begin
    n:=1;
  try
    for i:=1 to vle.RowCount-1 do
    begin
        n:=i;
        s:=vle.Cells[0,i];
        if  s='' then s:='0';
        Datax[i-1]:= strtofloat(s);
        s:=vle.Cells[1,i];
        if  s='' then s:='0';
        Datay[i-1]:= strtofloat(s);
    end;
  except
    application.MessageBox(Pchar('读取数据在读第'+inttostr(n)+'行数据时出现错误。'),'出错：',MB_ICONERROR);
  end;
end;

function  getDataMin(const Data:array of double;const n:integer):double;
var i       :integer;
    xTemp   :double;
begin
    xTemp:=Data[0];
    for i:=1 to n-1 do
    begin
        if xTemp>Data[i] then
            xTemp:=Data[i];
    end;  
    result:=xTemp;
end;

function  getDataMax(const Data:array of double;const n:integer):double;
var i       :integer;
    xTemp   :double;
begin
    xTemp:=Data[0];
    for i:=1 to n-1 do
    begin
        if xTemp<Data[i] then
            xTemp:=Data[i];
    end;
    result:=xTemp;
end;

function  getDataSum(const Data:array of double;const n:integer):double;
var i       :integer;
    xTemp   :double;
begin
    xTemp:=0;
    for i:=0 to n-1 do
    begin
        xTemp:=xTemp+Data[i];
    end;
    result:=xTemp;
end;

function  getDataProductSum(const Datax:array of double;const Datay:array of double;const n:integer):double;
var i       :integer;
    xTemp   :double;
begin
    xTemp:=0;
    for i:=0 to n-1 do
    begin
        xTemp:=xTemp+Datax[i]*Datay[i];
    end;
    result:=xTemp;
end;

function  getDataSr(const Datax:array of double;
                    const Datay:array of double;
                    const n:integer;
                    const A:double;
                    const B:double):double;
var i       :integer;
    xTemp   :double;
begin
    result:=0;
    xTemp:=0;
    try
        for i:=0 to n-1 do
        begin
            xTemp:=xTemp+sqr(Datay[i]-B*Datax[i]-A);
        end;
    except
        result:=Infinity;
        exit;
    end;

    if n<=2 then
    begin
        if xTemp=0 then
            result:=1
        else if xTemp>0 then
            result:=+Infinity
        else if xTemp<0 then
            result:=-Infinity;
    end
    else
    begin
        result:=sqrt(xTemp/(n-2));
    end;
end;

procedure Calculate(const Datax:array of double;const Datay:array of double;var xx:Txx);
var xTemp   :double;
    yTemp   :double;
begin
    xx.xMin:=getDataMin(Datax,xx.n);
    xx.yMin:=getDataMin(Datay,xx.n);
    xx.xMax:=getDataMax(Datax,xx.n);
    xx.yMax:=getDataMax(Datay,xx.n);
    xx.Hx:=getDataSum(Datax,xx.n);
    xx.Hy:=getDataSum(Datay,xx.n);
    xx.Hxy:=getDataProductSum(Datax,Datay,xx.n);
    xx.Hxx:=getDataProductSum(Datax,Datax,xx.n);
    xx.Hyy:=getDataProductSum(Datay,Datay,xx.n);
    xx.Px:=xx.Hx/xx.n;
    xx.Py:=xx.Hy/xx.n;
    // A
    xTemp:=xx.Hxx*xx.Hy-xx.Hx*xx.Hxy;
    yTemp:=xx.n*xx.Hxx-sqr(xx.Hx);
    if yTemp=0 then
    begin
        if xTemp=0 then
            xx.A:=Infinity
        else if xTemp>0 then
            xx.A:=+Infinity
        else if xTemp<0 then  
            xx.A:=-Infinity;
    end
    else
    begin
        xx.A:=xTemp/yTemp;
    end;

    // B
    xTemp:=xx.n*xx.Hxy-xx.Hx*xx.Hy;
    yTemp:=xx.n*xx.Hxx-sqr(xx.Hx);
    if yTemp=0 then
    begin
        if xTemp=0 then
            xx.B:=Infinity
        else if xTemp>0 then
            xx.B:=+Infinity
        else if xTemp<0 then  
            xx.B:=-Infinity;
    end
    else
    begin
        xx.B:=xTemp/yTemp;
    end;
    
    // r
    xTemp:=xx.Hxy-xx.Hx*xx.Hy/xx.n;
    yTemp:=sqrt((xx.Hxx-sqr(xx.Hx)/xx.n)*(xx.Hyy-sqr(xx.Hy)/xx.n));
    if yTemp=0 then
    begin
        if xTemp=0 then
            xx.r:=Infinity
        else if xTemp>0 then
            xx.r:=+Infinity
        else if xTemp<0 then
            xx.r:=-Infinity;
    end
    else
    begin
        xx.r:=xTemp/yTemp;
    end;
    xx.r:=abs(xx.r);//?


    //Sr
    xx.Sr:=getDataSr(Datax,Datay,xx.n,xx.A,xx.B);
    if abs(xx.Sr)=Infinity then
    begin
        xx.Sa:=Infinity;
        xx.Sb:=Infinity;
    end
    else
    begin
        // Sa
        yTemp:=xx.n*xx.Hxx-sqr(xx.Hx);
        if yTemp=0 then
            xx.Sa:=Infinity
        else
            xx.Sa:=sqrt(xx.Hxx/yTemp)*xx.Sr;
        // Sb
        yTemp:=xx.n*xx.Hxx-sqr(xx.Hx);
        if yTemp=0 then
            xx.Sb:=Infinity
        else
            xx.Sb:=sqrt(xx.n/yTemp)*xx.Sr;
    end;
end;

procedure DarwPic(var Img :TImage;const Datax:array of double;const Datay:array of double;const xx:Txx);
var X1,X2,Y1,Y2 :double;
    W,H         :integer;
    i           :integer;
    xt,yt       :integer;
    function getX(x:double):integer;
    begin
        result:=trunc((W-30)*(x-X1)/(X2-X1)+10);
    end;
    function getY(y:double):integer;
    begin
        result:=trunc((30-H)*(y-Y2)/(Y2-Y1)+10);
    end;
begin
    W:=img.Width;
    H:=img.Height;
    if xx.xMin=xx.xMax then
    begin
        X1:=xx.xMin-1;
        X2:=xx.xMax+1;
    end
    else
    begin
        X1:=xx.xMin;
        X2:=xx.xMax;
    end;
    if xx.yMin=xx.yMax then
    begin
        y1:=xx.yMin-1;
        y2:=xx.yMax+1;
    end
    else
    begin
        y1:=xx.yMin;
        y2:=xx.yMax;
    end;

    // 画线
    Img.Canvas.Pen.Color:=clLime; 
    img.Canvas.Pen.Width:=2;
    if abs(xx.B)=Infinity then
    begin
        img.Canvas.MoveTo(getX(xx.xMin),0);
        img.Canvas.LineTo(getX(xx.xMin),H);
    end
    else
    begin
        img.Canvas.MoveTo(getX(X1),getY(xx.A+xx.B*X1));
        img.Canvas.LineTo(getX(X2),getY(xx.A+xx.B*X2));
    end;
    
    // 画点
    Img.Canvas.Pen.Color:=clRed;
    img.Canvas.Pen.Width:=3;
    for i:=0 to xx.n-1 do
    begin
        xt:=getX(Datax[i]);
        yt:=getY(Datay[i]);
        img.Canvas.MoveTo(xt,yt);
        img.Canvas.LineTo(xt+1,yt);
    end;

end;

procedure SetVLEData(var VLE :TValueListEditor;const xx :Txx);overload;
begin
    Vle.Strings.Text:='';
    VleInsertRow(Vle,'数据组数n',xx.n);
    VleInsertRow(Vle,'线性回归方程系数A',xx.A);
    VleInsertRow(Vle,'线性回归方程系数B',xx.B);
    VleInsertRow(Vle,'相关系数r',xx.r);
    VleInsertRow(Vle,'标准误差Sr',xx.Sr);
    VleInsertRow(Vle,'A的标准误差Sa',xx.Sa);
    VleInsertRow(Vle,'B的标准误差Sb',xx.Sb);
    VleInsertRow(Vle,'x的平均值',xx.Px);
    VleInsertRow(Vle,'y的平均值',xx.Py);
    VleInsertRow(Vle,'x的和',xx.Hx);
    VleInsertRow(Vle,'y的和',xx.Hy);
    VleInsertRow(Vle,'x*y的和',xx.Hxy);
    VleInsertRow(Vle,'x*x的和',xx.Hxx);
    VleInsertRow(Vle,'y*y的和',xx.Hyy);
    VleInsertRow(Vle,'x的最小值xMin',xx.xMin);
    VleInsertRow(Vle,'x的最大值xMax',xx.xMax);
    VleInsertRow(Vle,'y的最小值yMin',xx.yMin);
    VleInsertRow(Vle,'y的最大值yMax',xx.yMax);
end;

procedure SetVLEData(var VLE :TValueListEditor);overload;
begin
    Vle.Strings.Text:='';
    Vle.InsertRow('数据组数n','',true);
    Vle.InsertRow('线性回归方程系数A','',true);
    Vle.InsertRow('线性回归方程系数B','',true);
    Vle.InsertRow('相关系数r','',true);
    Vle.InsertRow('标准误差Sr','',true);
    Vle.InsertRow('A的标准误差Sa','',true);
    Vle.InsertRow('B的标准误差Sb','',true);
    Vle.InsertRow('x的平均值','',true);
    Vle.InsertRow('y的平均值','',true);
    Vle.InsertRow('x的和','',true);
    Vle.InsertRow('y的和','',true);
    Vle.InsertRow('x*y的和','',true);
    Vle.InsertRow('x*x的和','',true);
    Vle.InsertRow('y*y的和','',true);
    Vle.InsertRow('x的最小值xMin','',true);
    Vle.InsertRow('x的最大值xMax','',true);
    Vle.InsertRow('y的最小值yMin','',true);
    Vle.InsertRow('y的最大值yMax','',true);
end;

procedure VLEInsertRow(var VLE:TValueListEditor;const str:string;const x:double);
var sTemp   :string;
begin
    if x=Infinity then
        sTemp:='无穷大'
    else if x=-Infinity  then
        sTemp:='负无穷大'
    else
        sTemp:=floattostr(x);
    Vle.InsertRow(str,sTemp,true);
end;

function  GetFileNamePath(const FileName:string):string;
var n   :integer;
    i   :integer;
    L   :integer;
    s   :string;
begin
    s:= FileName ;
    L:=length(s);
    n:=1;
    for i:=L downto 1 do
    begin
        n:=i;
        if s[i]='\' then break;
    end;
    if n=1 then
        s:='\'
    else
        setlength(s,n);
    result:=s;
end;

end.