unit CurveFitUnit;

interface

uses
  Windows, Messages, SysUtils, Variants, Classes, Graphics, Controls,
  Dialogs, Math;

type
	TDoubleXY = record
  	x: Double;
    y: Double;
  end;
  TDoubleArray = array of Double;

  procedure LSFit(const k: Integer; const XYList: array of TDoubleXY; var si,p,be,al,s: TDoubleArray; l: Boolean=true);

implementation

{ ***************************************************************************
Name: 最小二乘曲線擬合
Function: 本過程用于求列表函數(shù)在最小二乘意義下的最佳擬合廣義多項式和與此等價
      的冪多項式。它自動選擇擬合多項式的“最佳”次數(shù)。
ToDo: 設(shè)給定n個點 X1<X2<...<Xn 上的函數(shù)值 Y1,Y2,...,Yn 和與函數(shù)值的標(biāo)準(zhǔn)誤差平
      方成反比的權(quán) W1，W2，。。。，Wn，我們構(gòu)造 m（m≤n-1）次多項式:
              Ym(x) = C0P0(x) + c1P1(x) + ... + CmPm(x)                   (1)
      使
              (σm)^2 = ∑Wi(Ym(Xi)-Yi)^2,   i=1,2,...,n                  (2)
      取極小，這里 Pr(x) 是關(guān)于權(quán) Wi 滿足正交條件
              ∑WiPr(Xi)Ps(Xi) = 0,    r≠s, i=1,2,...,n                  (3)
      的 r 次多項式。

      導(dǎo)出正規(guī)方程組并利用正交條件(3)，得到:
              Cj = ∑WiYiPj(Xi) / ∑Wi(Pj(Xi))^2,   j=0,1,...,m           (4)
      正交多項式 Pj(x)可用如下遞推公式計算:
              P0(x) = 1,
              P1(x) = (x-α1)P0(x),
              Pj+1(x) = (x-αj+1)Pj(x) - βjPj-1(x),   j=1,2,...,m-1,
      其中
              αj+1 = (∑WiXi(Pj(Xi))^2)/(∑Wi(Pj(Xi))^2),
              βj   = (∑Wi(Pj(Xi))^2)/(∑Wi(Pj-1(Xi))^2).

      這樣，由(4)計算用正交多項式表示的廣義多項式的系數(shù)，再利用正交多項式的遞
      推公式計算 Pj(x) 之后，代入(1)并展開，就能得到與廣義多項式(1)等價的冪多
      項式:
              Q(x) = D0 + D1X + ... + DmX^m.

      結(jié)合光滑度和逼近度，尋找擬合多項式的“最佳”次數(shù) h(h≤m)。先計算平均平
      方誤差 δh = (σh)^2/(n-h-1)，如果當(dāng)擬合多項式次數(shù)逐次增加，h 是第一次
      使 δh≤δh+1 且 δj>0.6δh(m≥j>h+1) 的擬合多項式的次數(shù)，則 h 是“最
      佳”次數(shù)；如果 δh≤δh+1，但 δj≤0.6δh(j>h+1)，則“最佳”次數(shù)是j。

      另外，利用公式(1)，(2)，(3)和(4)可以導(dǎo)出計算 (σi)^2 的遞推公式:
              (σi)^2 = ∑Wj[Yi-1(Xj) + CiPi(Xj) - Yi]^2
                      = (σi-1)^2 - 2Ci∑WjPi(Xj)Yj + (Ci)^2 ∑Wj(Pi(Xj))^2
                      = (σi-1)^2 - (Ci)^2 ∑Wj(Pi(Xj))^2

      式中:
              j=1,2,...,n

Params:
      n       給定點的個數(shù)。
      k       表示用次數(shù)小于等于k的多項式進行擬合。
      x       數(shù)組x[0..n-1]，存放給定點，要求:
                  x[i] < x[i+1],    i=0,1,...,n-2。
      f       數(shù)組f[0..n-1]，存放給定點上的函數(shù)值。
      w       數(shù)組w[0..n-1]，存放給定點上的權(quán)。
      si      數(shù)組si[0..k]，存放平均平方誤差，即
                  si[j] = δj,    j=0,1,...,k。
      p       數(shù)據(jù)p[0..k]，存放冪多項式系數(shù)，p[i]是x^i的系數(shù)，i=0,1,..,k。
      l       布爾量。如果l為假，則自動選擇“最佳”次數(shù)；如果l為真，則用k次
              多項式進行擬合。
      al      數(shù)組al[0..k-1]，存放αi。
      be      數(shù)組be[0..k]，存放βi。
      s       數(shù)組s[0..k]，存放廣義多項式的系數(shù)Ci。

      說明:   為了方便使用，這里采用了參數(shù) XYList 來代替 x和f 數(shù)組存放給定點
              和給定點上的函數(shù)值，給定點的個數(shù) n 已經(jīng)隱含在參數(shù) XYList 中了。
              給定點上的權(quán)數(shù)組 w 恒設(shè)為1。
              調(diào)用本函數(shù)時，用戶只需聲明 si,p,be,al,s 參數(shù)，可以不為其分配存
              儲空間，本函數(shù)內(nèi)部將自動為它們分配適當(dāng)?shù)拇鎯臻g。

              請務(wù)必保證---- n ≥k+2

*************************************************************************** }

procedure PolyX(a,b: Double; var c,d: TDoubleArray; m,n: Integer);
var
  i,j,k: Integer;
  z,w: TDoubleArray;
begin
  SetLength(z,m+1);
  SetLength(w,m+1);

  z[0] := 1.0;
  w[0] := 1.0;
  d[0] := c[0];

  for i := 1 to n do
  begin
    w[i] := 1.0;
    z[i] := b * z[i-1];
    d[0] := d[0] + c[i] * z[i];
  end;

  for j := 1 to n do
  begin
    w[0] := w[0] * a;
    d[j] := c[j] * w[0];
    k := 1;
    for i := j+1 to n do
    begin
      w[k] := a * w[k] + w[k-1];
      d[j] := d[j] + c[i] * w[k] * z[k];
      Inc(k);
    end;
  end;
end;

procedure LSFit(const k: Integer; const XYList: array of TDoubleXY; var si,p,be,al,s: TDoubleArray; l: Boolean=true);
var
  i,j,n: Integer;
  w: TDoubleArray;
  ctp,cpsave,cp: TDoubleArray;
  lp,tp: TDoubleArray;
  clp: TDoubleArray;
  du,delsq,om,lw,tw,simin,a,b: Double;
  swx,comp: Boolean;

  x,f: TDoubleArray;
Label
  L1;
begin
  n := Length(XYList);
  SetLength(x,n);             //0~n-1
  SetLength(f,n);
  SetLength(w,n);
  for i := 0 to n-1 do
  begin
    x[i] := XYList[i].x;
    f[i] := XYList[i].y;
    w[i] := 1.0;
  end;

  SetLength(ctp,k+1);         //0~k
  SetLength(cpsave,k+1);
  SetLength(cp,k+1);

  SetLength(lp,n);            //0~n-1
  SetLength(tp,n);

  SetLength(clp,k+2);

  SetLength(si,k+1);          //0~k
  SetLength(p,k+1);
  SetLength(be,k+1);
  SetLength(al,k);            //0~k-1
  SetLength(s,k+1);           //0~k

  for i := 0 to k do cp[i] := 0;
  simin := 0;
  swx := true;
  be[0] := 0;
  clp[0] := 0;
  clp[1] := 0;                //-1
  delsq := 0;
  om := 0;
  ctp[0] := 1.0;
  tw := 0;
  comp := false;

  for i := 0 to n-1 do        //1-n
  begin
    tp[i] := 1.0;
    lp[i] := 0;
    delsq := delsq + w[i] * Power(f[i],2);
    om := om + w[i] * f[i];
    tw := tw + w[i];
  end;

  cp[0] := om / tw;
  s[0] := cp[0];
  delsq := delsq - s[0] * om;
  si[0] := delsq / (n-1);
  a := 4.0 / (x[n-1]-x[0]);   //n,1
  b := -2 - a * x[0];

  for i := 0 to n-1 do        //1-n
    x[i] := a * x[i] + b;

  for i := 0 to k-1 do
  begin
    du := 0;
    for j := 0 to n-1 do      //1-n
      du := du + w[j] * x[j] * Power(tp[j],2);
    al[i] := du / tw;
    lw := tw;
    tw := 0;
    om := 0;

    for j := 0 to n-1 do      //1-n
    begin
      du := be[i] * lp[j];
      lp[j] := tp[j];
      tp[j] := (x[j] - al[i]) * tp[j] - du;
      tw := tw + w[j] * Power(tp[j],2);
      om := om + w[j] * f[j] * tp[j];
    end;

    be[i+1] := tw / lw;
    s[i+1] := om / tw;
    delsq := delsq - s[i+1] * om;
    si[i+1] := delsq / (n-i-1);

    if l then goto L1;
    if not comp then
    begin
      if swx then
      begin
        if si[i+1] >= si[i] then
        begin
          swx := false;
          simin := si[i];
          for j := 0 to k do
            cpsave[j] := cp[j];
        end;
        goto L1;
      end;

      if si[i+1] < 0.6*simin then comp := true;
L1:
      for j := 0 to i do
      begin
        du := clp[j] * be[i];
        clp[j] := ctp[j];
        if j = 0 then
          ctp[j] := 0 - al[i] * ctp[j] - du
        else
          ctp[j] := clp[j-1] - al[i] * ctp[j] - du;
        cp[j] := cp[j] + s[i+1] * ctp[j];
      end;

      cp[i+1] := s[i+1];
      ctp[i+1] := 1.0;
      clp[i+1] := 0;
      if (not (comp or swx)) and (i = k-1) then
      begin
        for j := 0 to k do
          cp[j] := cpsave[j];
      end;
    end;
  end;

  PolyX(a,b,cp,p,n,k);
end;

end.