データの分布を正規分布に近づけるための変換を行います。１次元の最適化
計算を行う必要があり、最適化ルーチンには黄金分割法を使いました。
　Box-Cox変換の詳細は統計学の参考書を、黄金分割法については数値計算の
参考書を参照してください。

{================== 以下がソース ======================}
uses SysUtils, Math;

Type
  (*----------------------------------------------
    不定サイズの配列を動的に確保するための型
    動的配列が利用できるときにはこんなものは不要
  -----------------------------------------------*)
  TVector = array[0..0] of extended;
  PVector = ^TVector;

(*--------------------------------------------------------------------------
  データのBox-Cox変換
  与えられたデータの配列を変換して返す。
  関数の値としては変換時のλの値を返す

  Box-Cox変換の内容

  Y' = (Y^λ - 1) / λ (λ<>0)
  Y' = lnY             (λ= 0)
  L = -ν/2 ln(st2) + (λ-1)*ν/n*�罵nY

  Y  : データ
  Y' : 変換したデータ
  ν : 自由度
  n  : データ数
  st2: 変換したデータの分散（標本分散）

  尤度Lを最大化するλの値で変換すると、正規分布に近づけるのに最もよい
  変換になる。λの最適化ルーチンには黄金分割法を使用。

  参考文献
　Box-Cox変換については、
  Sokal & Rohlf(1995) "Biometry 3rd Ed." pp.417-418.
  黄金分割法による最適化については
  Shammas(1995)「C/C++数値計算アルゴリズムブック」
  　　鈴木治郎訳(1997) 海文堂 pp.123-124.
-------------------------------------------------------------------------*)
function BoxCoxTransformation(var X : array of extended) : extended;
var
  lambda : extended;  (* λ                   *)
  SlnY   : extended;  (* �罵nY                *)
  n,DF   : integer;   (* n, ν                *)
  Y,Yt   : PVector;   (* 配列確保用のポインタ *)
  i      : integer;   (* カウンタ             *)

  (* 以下は黄金分割用の変数 *)
  a,b,c,d,
  Fc,Fd  : extended;
  t      : extended;       (* 黄金分割比     *)
  Tol    : extended;       (* 許容誤差       *)
  iteration    : integer;  (* 反復回数       *)
  MaxIteration : integer;  (* 反復回数の上限 *)

  {最小化すべき尤度関数ただし符号を逆にしてある}
  function L(r : extended) : extended;
  var
    i                   : integer;
    SX,Average,Variance : extended;
  begin
    (* 変換データの計算 *)
    if r <> 0 then
      for i := Low(X) to High(X) do Yt^[i] := (Power(Y^[i], r) - 1) / r
    else
      for i := Low(X) to High(X) do Yt[i] := ln(Y[i]);
    (* Ytの標本分散の計算 *)
    Sx := 0; Variance := 0;
    for i := Low(X) to High(X) do SX  := SX + Yt^[i];
    Average := SX / n;
    for i := Low(X) to High(X) do
      Variance := Variance + Sqr(Average - Yt^[i]);
    Variance := Variance / DF;
    (* 尤度の計算(符号反転) *)
    Result := DF/2 * ln(Variance) - (r - 1) * DF / n * SlnY;
  end;

begin
  n := High(X) + 1;
  DF:= High(X);
  (* ワーク用の配列メモリの確保 *)
  GetMem(Y , SizeOf(extended) * n);
  GetMem(Yt, SizeOf(extended) * n);
  for i := Low(X) to High(X) do Y^[i] := X[i];
  SlnY := 0;
  for i := Low(X) to High(X) do SlnY := SlnY + ln(X[i]);

  (* 以下はλの最適化ルーチン(黄金分割法)        *)
  Tol := 1.0E-10;          (* 許容誤差の設定     *)
  MaxIteration := 100;     (* 最大反復回数の設定 *)
  t := (Sqrt(5)-1)/2;      (* 黄金分割比         *)
  iteration := 0;
  a := -100;
  b := 100;
  c := a + (1 - t) * (b - a);
  d := b - (1 - t) * (b - a);
  Fc := L(c);
  Fd := L(d);
  while (Abs(b - a) > Tol) and (Iteration < MaxIteration) do begin
    Inc(Iteration);
    if Fc < Fd then begin
      b := d;
      d := c;
      c := a + (1 - t) * (b - a);
      Fd := Fc;
      Fc := L(c);
    end else begin
      a := c;
      c := d;
      d := b - (1 - t) * (b - a);
      Fc := Fd;
      Fd := L(d);
    end;
  end;
  if Iteration < MaxIteration then begin
    lambda := (a + b) / 2;
    if lambda <> 0 then
      for i := Low(X) to High(X) do
        X[i] := (Power(Y^[i], lambda) - 1) / lambda
    else
      for i := Low(X) to High(X) do
        X[i] := ln(Y[i]);
    Result := lambda;
  end else
    Raise Exception.Create('too many iterations in BoxCoxTransformation');
  (* メモリの開放 *)
  FreeMem(Y);
  FreeMem(Yt);
end;

                              2000/12/19(Tue) 04:46pm  GFD03044 ひの
 


- FDELPHI  MES(16):玉石混淆みんなで作るSample蔵【見本蓄積】 00/12/27 -

Original document by ひの            氏 ID:(GFD03044)

ここにあるドキュメントは NIFTY SERVEの Delphi Users' Forum の16番会議室「玉石混淆みんなで作るSample蔵」に投稿されたサンプルです。これらのサンプルはボーランド株式会社がサポートする公式のものではありません。また、必ずしも動作が検証されているものではありません。これらのサンプルを使用したことに起因するいかなる損害も投稿者、およびフォーラムスタッフはその責めを負いません。使用者のリスクの範疇でご使用下さい。