Statistics::OLS - perform ordinary least squares and associated statistics, v 0.07.

Location : Home > Languages > Perl > Package
Title : Statistics::OLS

　Statistics::OLS - 最小二乗法

use Statistics::OLS; 

my $ls = Statistics::OLS->new(); 

$ls->setData (\@xydataset) or die( $ls->error() );
$ls->setData (\@xdataset, \@ydataset);

$ls->regress();

my ($intercept, $slope) = $ls->coefficients();
my $R_squared = $ls->rsq();
my ($tstat_intercept, $tstat_slope) = $ls->tstats();
my $sigma = $ls->sigma();
my $durbin_watson = $ls->dw();  

my $sample_size = $ls->size();
my ($avX, $avY) = $ls->av();
my ($varX, $varY, $covXY) = $ls->var();
my ($xmin, $xmax, $ymin, $ymax) = $ls->minMax();

# setData() で最初に呼び出された値に依存する x-y データまた y データを返す。
my @predictedYs = $ls->predicted();
my @residuals = $ls->residuals();

説明

　私は２次元データ y = a + bx における最小二乗法（線形曲線フィティング）を実装するために Statistics::OLS を書いた。
　CPAN で見つけた統計モジュール（Statistics::Descriptive）は１変数の解析を目的に作成されてたものである。これには OLS も添付されているが、幾分柔軟性に欠け、二変数解析が不十分であった。しかしいつかNevertheless, it might make sense to fold OLS をそのモジュールに組み込むか、上位互換のものを作成することは意味のあることである。
　Statistics::OLS は回帰直線の傾きと切片、 T 統計値、決定係数（R squared）、回帰の標準誤差、ダービン・ワトソン統計量（Durbin-Watson statistic）を計算する。また残差も返す。
　２次元の最小２乗法の処理は非常にシンプルであるが、重回帰は非常に難しい。そのため OLS は複数の独立変数に対しては十分に機能しない。
　これはベータコードであり、十分にテストされていない。いくつかの公開されているデータセットでは確認している。フィードバックは歓迎する。特に、エラーや正しく計算されない結果を教えてくれることはありがたい。

使用法

新たな regression オブジェクトを生成：new()

use Statistics::OLS;
my $ls = Statistics::OLS->new;

データセットを登録：setData()

$ls->setData (\@xydata);
$ls->setData (\@xdata, \@ydata);

　setData() メソッドは２次元のデータセットを regression オブジェクトに登録する。データセットは x,y データの対を含む配列への参照としても、xデータ・yデータの２つの配列への参照としても設定することができる。（どちらの場合もスクリプト内のデータ配列はキャッシュ－オブジェクト内へのコピー－はされない。）　もしデータを変更したら再度 setData() をオプションで呼び出すことができるが（エラーをチェックしたければ下記を確認のこと）、少なくとも regress() メソッド（下記参照）は新しいデータに基づく再計算のために呼び出さなければならない。または単純にデータを変更しないこと。
　スクリプト内で１次元配列の場合、 n+1 個の x,y データ点を含む (x0, y0, ..., xn, yn) という形式の配列をコンストラクトする。それから setData() メソッドにデータ配列への参照を渡す。（もし配列とは何かを知らなければ、 setData() への引数として渡すときにデータ配列の名前の前にバックスラッシュ（\）を置くだけでよい。
　このような感じだ。

my @xydata = qw( -3 9   -2 4   -1 1   0 0   1 1  2 4  3 9);
$ls->setData (\@xydata);

　または、一方は水平方向、他方は対応する垂直方向のデータを格納した、２つの同じ長さの配列をコンストラクトするほうが便利かも知れない。それから両方の配列（ただし水平方向のデータを先）を setData() に渡す。

my @xdata = qw( -3  -2  -1  0  1  2  3 );
my @ydata = qw(  9   4   1  0  1  4  9 );
$ls->setData (\@xdata, \@ydata);

　エラーチェック： setData() メソッドは成功時には正整数を返し、失敗時には 0 を返す。 setData() が落ちた場合には error() メソッドから出力されるエラーメッセージを取得することができる。エラーメッセージの文字列は以下のうちの１つである。

The data set does not contain an equal number of x and y values. 
The data element ... is non-numeric.
The data set must contain at least three points.

　スクリプト内ではエラーに関して以下のようにテストすればよい。

$ls->setData (\@data) or die( $ls->error() );

　現在のバージョンでは数値と小数点（ヨーロッパの皆さん、ごめんなさい）、科学的表記（1.7E10）、負の符号のみに対応している。通貨記号（$）、時刻（11:23am）や日付（23/05/98）はサポートしていないし、エラーが出力される。これはいずれ対応するかも知れない。

回帰分析：regress()

$ls->regress() or die ( $ls->error() );

　これが計算のほとんどである。データを設定してからこのメソッドを呼び出す。しかし回帰結果を呼び出す前に、である。データを変更すれば過去の計算は不正確になるため、再度呼び出さなければならない。regress() メソッドは成功時には 1 を返す。エラーメッセージは

No datset has been registered.

のみである。
　後述するような特別な統計量に関しては（ゼロ割などの影響で） undef を返すこともある。

回帰結果の取得：coefficients(), rsq(), tstats() 等

my ($intercept, $slope) = $ls->coefficients();
my $R_squared = $ls->rsq();
my ($tstat_intercept, $tstat_slope) = $ls->tstats();
my $sigma = $ls->sigma();
my $durbin_watson = $ls->dw();

my $sample_size = $ls->size();
my ($avX, $avY) = $ls->av();
my ($varX, $varY, $covXY) = $ls->var();
my ($xmin, $xmax, $ymin, $ymax) = $ls->minMax();

　regress() を呼び出した後にこれらのメソッドを呼び出す。これらのほとんどは計量経済学ではなじみ深いものである。傾きが無限大（X の変動幅が 0 ）であれば undef となる。決定係数（R-squared）は X または Y の標本分散が０（またはデータが共線形である）場合に 1.0 となる。)　もし X の分散が 0 であれば、T 統計量も undef となる。sigma は誤差項の等分散標準偏差推定値で推定標準偏差として知られる。 n-1 で計算している。Durbin-Watson はデータが共線形であれば undef を返す。

予測値または残差を返す：predicted(), residuals()

my @predictedYs = $ls->predicted();
my @residuals = $ls->residuals();

　regress() を呼び出した後でのみこれらのメソッドを呼び出す。双方のメソッドとも setData() で用いたものと同じフォーマットの配列を返す。データが setData() に (x0, y0, ..., xn, yn) という形式の @xydata 配列への参照を渡されたらこれらのメソッドの結果は同じ形式で返される。ただし y 値が係数から推計された予測値であるか, or the 観測値と予測値の間の残差誤差である場合を除く。
　データが２つの配列 @xdata = (x0 ... xn) 及び @ydata = (y0 ... yn) への参照として渡された場合、これら２つのメソッドの結果は y タイプデータの１つの配列となる。予測値か残差誤差のいずれかである。元の x データ配列はこれらの結果配列と対応している。

バグと To Do

　本モジュールはベータコードであり、正しく稼動することを保証しない。私はテストを十分にこなせていない。
　将来的には、日付・時刻・通貨など対応データフォーマットの拡張を行う。
　重回帰への拡張はおそらくありえない。非常に難しいからだ。Fortran や Perl Data Language に基づいたほかのものを使ったほうがいいだろう。
　これをさらに包括的なライブラリとして Statistics::Descriptive に組み込むことは意味があるだろう。しかしそれは近未来のことではない。大プロジェクトになりそうだからだ。
　コメント及びバグ報告を歓迎する。

著者

　本プログラムはフリーソフトウェアであり、Perl 本体と同等の条件で修正／再配布してもよい。
　この作業は、要望を出してくれた故 Dr. Andrew Morton に捧げられる。安らかに眠れ、友よ。

参考資料

　Statistics::Descriptive(1) モジュールは有用な一変数統計を扱っており、十分にテストされている。Perl Data Language （CPAN を見よ）も統計を扱うには十分であることが確認されている。
　線形回帰についてはあらゆる計量経済学やほとんどの確率・統計のテキストで議論されている。私は Basic Econometrics 2nd ed., by Gujaratii, New York: McGraw-Hill,1988, を用いており、公式やテストサンプルにはここからとった。（appendix 3A.6, page 87）.

Updated : 2007/05/15