unitm2.pas

Delphi编写的最小二乘法的曲线拟合算法

unit UnitM2;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls, Forms,
  Dialogs, StdCtrls, ExtCtrls, Spin;

type
  TForm1 = class(TForm)
    Memo1: TMemo;
    Splitter1: TSplitter;
    Button1: TButton;
    Panel1: TPanel;
    Panel2: TPanel;
    Image1: TImage;
    Label1: TLabel;
    SpinEdit1: TSpinEdit;
    Label2: TLabel;
    Label3: TLabel;
    Edit1: TEdit;
    Edit2: TEdit;
    Edit3: TEdit;
    procedure Button1Click(Sender: TObject);
  private

  public

  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses
  Math;

type
  TFloatAy=array of Double;
{
 //X,Y --  X,Y两轴的坐标
 //M   --  结果变量组数
 //N   --  采样数目
 //A   --  结果参数
}
function CalculateCurveParameter(X, Y: TFloatAy; M, N: Integer; A: TFloatAy):Boolean;
var
  i,j,k:Integer;
  Z,D1,D2,C,P,G,Q: Double;
  B,T,S:TFloatAy;
begin
  Result:=false;
  if N<1 then exit;
  SetLength(B,N);
  SetLength(T,N);
  SetLength(S,N);
  if M>N then
    M:=N;
  for i:=0 to Pred(M) do
    A[i]:=0;
  Z:=0;
  B[0]:=1;
  D1:=N;
  P:=0;
  C:=0;
  for i:=0 to Pred(N) do
  begin
    P:=P+X[i]-Z;
    C:=C+Y[i];
  end;
  C:=C/D1;
  P:=P/D1;
  A[0]:=C*B[0];
  if M>1 then
  begin
    T[1]:=1;
    T[0]:=-P;
    D2:=0;
    C:=0;
    G:=0;
    for i:=0 to Pred(N) do
    begin
      Q:=X[i]-Z-P;
      D2:=D2+Q*Q;
      C:=Y[i]*Q+C;
      G:=(X[i]-Z)*Q*Q+G;
    end;
    C:=C/D2;
    P:=G/D2;
    Q:=D2/D1;
    D1:=D2;
    A[1]:=C*T[1];
    A[0]:=C*T[0]+A[0];
    for j:=2 to Pred(M) do
    begin
      S[j]:=T[j-1];
      S[j-1]:=-P*T[j-1]+T[j-2];
      if j>=3 then
      begin
        for k:=j-2 downto 1 do
          S[k]:=-P*T[k]+T[k-1]-Q*B[k];
      end;
      S[0]:=-P*T[0]-Q*B[0];
      D2:=0;
      C:=0;
      G:=0;
      for i:=0 to Pred(N) do
      begin
        Q:=S[j];
        for k:=j-1 downto 0 do
          Q:=Q*(X[i]-Z)+S[k];
        D2:=D2+Q*Q;
        C:=Y[i]*Q+C;
        G:=(X[i]-Z)*Q*Q+G;
      end;
      C:=C/D2;
      P:=G/D2;
      Q:=D2/D1;
      D1:=D2;
      A[j]:=C*S[j];
      T[j]:=S[j];
      for k:=j-1 downto 0 do
      begin
        A[k]:=C*S[k]+A[k];
        B[k]:=T[k];
        T[k]:=S[k];
      end;
    end;
  end;
  Result:=true;
end;

procedure TForm1.Button1Click(Sender: TObject);
var
  x,y,a:TFloatAy;
  i,j,m,n:Integer;
  SL:TStringList;
  MinX,MaxX,MinY,MaxY:Double;
  px,py,mw,mh:Integer;
  fx,fy:Double;
begin
  Memo1.Text:=Trim(Memo1.Text);
  n:=Memo1.Lines.Count;
  if n<2 then exit;
  if Memo1.Lines[n-1]='' then
    Dec(n);
  if n<2 then exit;
  m:=SpinEdit1.Value;
  if m<1 then
    m:=1;
  if n<m then
  begin
    MessageDlg('点数不够！', mtError, [mbOK], 0);
    exit;
  end;
  SetLength(x,n);
  SetLength(y,n);
  SetLength(a,m);
  SL:=TStringList.Create;
  for i:=0 to Pred(n) do
  begin
    SL.CommaText:=Memo1.Lines[i];
    if SL.Count>1 then
    begin
      x[i]:=StrToFloatDef(SL[0],0);
      y[i]:=StrToFloatDef(SL[1],0);
    end
    else begin
      y[i]:=StrToFloatDef(Memo1.Lines[i],0);
      x[i]:=StrToFloatDef(Edit1.Text,0)+i*StrToFloatDef(Edit2.Text,1);
    end;
    if i=0 then
    begin
      MinX:=x[i];
      MaxX:=x[i];
      MinY:=y[i];
      MaxY:=y[i];
    end
    else begin
      if x[i]<MinX then MinX:=x[i];
      if x[i]>MaxX then MaxX:=x[i];
      if y[i]<MinY then MinY:=y[i];
      if y[i]>MaxY then MaxY:=y[i];
    end;
  end;
  SL.Free;
  if MinX=MaxX then exit;
  if MinY=MaxY then
  begin
    MinY:=MinY-1;
    MaxY:=MaxY+1;
  end;
  fy:=MaxY-MinY;
  MinY:=MinY-0.1*fy;
  MaxY:=MaxY+0.1*fy;
  //绘图
  with Image1 do
  begin
    Brush.Color:=clWhite;
    Brush.Style:=bsSolid;
    Canvas.FillRect(Rect(0,0,Width,Height));
  end;
  mw:=Image1.Width-2;
  mh:=Image1.Height-2;
  for i:=0 to Pred(n) do
  begin
    px:=Round(mw*(x[i]-MinX)/(MaxX-MinX))+1;
    py:=mh-Round(mh*(y[i]-MinY)/(MaxY-MinY))+1;
    Image1.Canvas.Pixels[px,py]:=clRed;
    Image1.Canvas.Pixels[px-1,py]:=clMaroon;
    Image1.Canvas.Pixels[px+1,py]:=clMaroon;
    Image1.Canvas.Pixels[px,py+1]:=clMaroon;
    Image1.Canvas.Pixels[px,py-1]:=clMaroon;
  end;
  Application.ProcessMessages;
  CalculateCurveParameter(x,y,m,n,a);
  Edit3.Text:=Format('y = %f',[a[0]]);
  for i:=1 to Pred(m) do
    Edit3.Text:=Edit3.Text+Format(' + %.8e * x^%d',[a[i],i]);
  Application.ProcessMessages;
  for i:=0 to Pred(mw) do
  begin
    fx:=(MaxX-MinX)*(i/mw)+MinX;
    fy:=0;
    for j:=Pred(m) downto 1 do
      fy:=fy+Power(fx,j)*a[j];
    fy:=fy+a[0];
    px:=i+1;
    py:=mh-Round(mh*(fy-MinY)/(MaxY-MinY))+1;
    Image1.Canvas.Pixels[px,py]:=clGreen;
  end;
  SetLength(x,0);
  SetLength(y,0);
  SetLength(a,0);
end;

end.