最小二乗法でy=ax+bの回帰直線を求める ~Delphiでお手軽プログラミング

最小二乗法でy=ax+bの回帰直線を求める ~Delphiでお手軽プログラミング

最小二乗法のユニット

ファイル名「ULeastSquares.pas」で保存してください。
unit ULeastSquares;
{**************************************************************
  最小2乗法 ユニット
  y = ax + b
  で、x,yの配列(TLeastSquaresValues)を与えると、
  傾きTLeastSquaresAns.aと切片TLeastSquaresAns.bを返す
**************************************************************}

interface

type
  TLeastSquaresValue=record
    x,y:Extended;
  end;
  TLeastSquaresValues=TArray<TLeastSquaresValue>;
  TLeastSquaresAns=record
    a,b:Extended;
  end;
  TLeastSquares=class
    public
      class function Calc(XYArr:TLeastSquaresValues):TLeastSquaresAns;
    private
      class function SumXY(XYArr:TLeastSquaresValues):Extended;
      class function SumX(XYArr:TLeastSquaresValues):Extended;
      class function SumY(XYArr:TLeastSquaresValues):Extended;
      class function SumX2(XYArr:TLeastSquaresValues):Extended;
  end;

implementation

class function TLeastSquares.Calc(XYArr: TLeastSquaresValues): TLeastSquaresAns;
var x,y,xy,x2:Extended;
    e:Extended;
    n:Extended;
begin
  n:=Length(XYArr);//xyデータの数
  x:=SumX(XYArr);  //Σx
  y:=SumY(XYArr);  //Σy
  xy:=SumXY(XYArr);//Σx*y
  x2:=SumX2(XYArr);//Σx*x
  e:=(n*x2-x*x);
  if e<>0 then
  begin
    result.a:=(n*xy-x*y)/e;
    result.b:=(x2*y-xy*x)/e;
  end
  else
  begin
    result.a:=0;
    result.b:=0;
  end;
end;

class function TLeastSquares.SumX(XYArr: TLeastSquaresValues): Extended;
var i:Integer;
begin
  //Σx
  result:=0;
  for i := Low(XYArr) to High(XYArr) do
    result:=result+XYArr[i].x;
end;

class function TLeastSquares.SumY(XYArr: TLeastSquaresValues): Extended;
var i:Integer;
begin
  //Σy
  result:=0;
  for i := Low(XYArr) to High(XYArr) do
    result:=result+XYArr[i].y;
end;

class function TLeastSquares.SumXY(XYArr: TLeastSquaresValues): Extended;
var i:Integer;
begin
  //Σx*y
  result:=0;
  for i := Low(XYArr) to High(XYArr) do
    result:=result+XYArr[i].x*XYArr[i].y;
end;

class function TLeastSquares.SumX2(XYArr: TLeastSquaresValues): Extended;
var i:Integer;
begin
  //Σx*x
  result:=0;
  for i := Low(XYArr) to High(XYArr) do
    result:=result+XYArr[i].x*XYArr[i].x;
end;

end.

使用例

ファイル⇒新規作成⇒VCLフォームアプリケーションをクリックします(新規プロジェクトが作成されます)。
フォームにTButton×1個、TMemo×1個をドラッグ&ドロップします。

Button1をダウブルクリックして以下ソースコードを記述します。
unit Unit1;

interface

uses
  Winapi.Windows, Winapi.Messages, System.SysUtils, System.Variants, System.Classes,
  Vcl.Graphics, Vcl.Controls, Vcl.Forms, Vcl.Dialogs, Vcl.StdCtrls;

type
  TForm1 = class(TForm)
    Button1: TButton;
    Memo1: TMemo;
    procedure Button1Click(Sender: TObject);
  private
    { Private 宣言 }
  public
    { Public 宣言 }
  end;

var
  Form1: TForm1;

implementation

{$R *.dfm}

uses ULeastSquares;

procedure TForm1.Button1Click(Sender: TObject);
var xy:TLeastSquaresValues;
    ans:TLeastSquaresAns;
begin
  SetLength(xy,10);
  xy[0].x:= 0.6094;	xy[0].y:= 0.8393;
  xy[1].x:=-0.2344;	xy[1].y:= 0.0354;
  xy[2].x:=-0.6029;	xy[2].y:=-0.1950;
  xy[3].x:= 0.1168;	xy[3].y:= 0.2681;
  xy[4].x:=-0.7346;	xy[4].y:=-0.8405;
  xy[5].x:=-0.8232;	xy[5].y:= 0.4612;
  xy[6].x:= 0.6606;	xy[6].y:=-0.3923;
  xy[7].x:= 0.926 ;	xy[7].y:= 0.9928;
  xy[8].x:=-0.5316;	xy[8].y:= 0.5900;
  xy[9].x:=-0.8263;	xy[9].y:=-0.0484;

  ans:=TLeastSquares.Calc(xy);

  Memo1.Lines.Add(
    Format('y = %8.4fx+%8.4f',[ans.a,ans.b])
  );
end;

end.

実行する

実行して、Button1を押すと、最小二乗法で傾きと切片が計算されます。

参考:エクセルでの近似直線結果