curvefitunit.pas

最小二乘相关介绍,最小二乘相关介绍,及其原程序!仔细整理收获不少!

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: 本过程用于求列表函数在最小二乘意义下的最佳拟合广义多项式和与此等价
      的幂多项式。它自动选择拟合多项式的“最佳”次数。
ToDo: 设给定n个点 X1<X2<...<Xn 上的函数值 Y1,Y2,...,Yn 和与函数值的标准误差平
      方成反比的权 W1，W2，。。。，Wn，我们构造 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) 是关于权 Wi 满足正交条件
              ∑WiPr(Xi)Ps(Xi) = 0,    r≠s, i=1,2,...,n                  (3)
      的 r 次多项式。

      导出正规方程组并利用正交条件(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)计算用正交多项式表示的广义多项式的系数，再利用正交多项式的递
      推公式计算 Pj(x) 之后，代入(1)并展开，就能得到与广义多项式(1)等价的幂多
      项式:
              Q(x) = D0 + D1X + ... + DmX^m.

      结合光滑度和逼近度，寻找拟合多项式的“最佳”次数 h(h≤m)。先计算平均平
      方误差 δh = (σh)^2/(n-h-1)，如果当拟合多项式次数逐次增加，h 是第一次
      使 δh≤δh+1 且 δj>0.6δh(m≥j>h+1) 的拟合多项式的次数，则 h 是“最
      佳”次数；如果 δh≤δh+1，但 δj≤0.6δh(j>h+1)，则“最佳”次数是j。

      另外，利用公式(1)，(2)，(3)和(4)可以导出计算 (σ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       给定点的个数。
      k       表示用次数小于等于k的多项式进行拟合。
      x       数组x[0..n-1]，存放给定点，要求:
                  x[i] < x[i+1],    i=0,1,...,n-2。
      f       数组f[0..n-1]，存放给定点上的函数值。
      w       数组w[0..n-1]，存放给定点上的权。
      si      数组si[0..k]，存放平均平方误差，即
                  si[j] = δj,    j=0,1,...,k。
      p       数据p[0..k]，存放幂多项式系数，p[i]是x^i的系数，i=0,1,..,k。
      l       布尔量。如果l为假，则自动选择“最佳”次数；如果l为真，则用k次
              多项式进行拟合。
      al      数组al[0..k-1]，存放αi。
      be      数组be[0..k]，存放βi。
      s       数组s[0..k]，存放广义多项式的系数Ci。

      说明:   为了方便使用，这里采用了参数 XYList 来代替 x和f 数组存放给定点
              和给定点上的函数值，给定点的个数 n 已经隐含在参数 XYList 中了。
              给定点上的权数组 w 恒设为1。
              调用本函数时，用户只需声明 si,p,be,al,s 参数，可以不为其分配存
              储空间，本函数内部将自动为它们分配适当的存储空间。

              请务必保证---- 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.