Title : Statistics::GaussHelmert
|
Location : Home > Languages > Perl > Package Title : Statistics::GaussHelmert |
|
Statistics::GaussHelmert - 汎用重み付け最小二乗法による推計
本ドキュメントは version 0.03 である。
use Statistics::GaussHelmert;
# 空のモデルを生成
my $estimation = new Statistics::GaussHelmert;
# 所与のモデル観測値 $y, 共分散行列を設定
# $Sigma_yy, an initial guess $b0 for the unknown parameters.
$estimation->observations($y);
$estimation->covariance_observations($Sigma_yy);
$estimation->initial_guess($b0);
# 陰モデル関数その閉包を使ってヤコビアン行列を指定
$estimation->observation_equations(sub { ... });
$estimation->Jacobian_unknowns(sub { ... });
$estimation->Jacobian_observations(sub { ... });
# 必須ではないが未知のパラメータに対する条件を付与したければ
$estimation->constraints(sub { ... });
$estimation->Jacobian_constraints(sub { ... });
# 推定開始
$estimation->start(verbose => 1);
# 結果出力
print $estimation->estimated_unknown(),
$estimation->covariance_unknown();
本モジュールは所与の観測値に対しモデルのパラメータを推定するための柔軟なツールである。モジュールは線形の推定モデルを提供し、基礎としているのはいわゆるガウス＝ヘルマートモデル（Gauss-Helmert model）である。
Statistics::GaussHelmert は Statistics::OLS とは任意の次元の任意の関数に適合させる点において異なっている。最小化する関数を陰関数で与えなければならない（伝統的な回帰モデルでは用関数で指定する）。未知数及び観測値に関しては微分を与える。未知数に関する条件の関数を微分で与える。さらにある程度適切な解が必要である。いくつかの問題では、いくつかの観測値と未知数で直接解くよりも適切な解を求めることが容易なことがある。
モデルの設定
知りたい未知のパラメータに関連したいくつかの観測値があるとしよう。たとえば２次元の観測値が未知の直線上状にあるものとする。
２次元上の観測値は PDL::Matrix オブジェクト $y （いわゆる multipiddle すなわち multi-PDL オブジェクト、すなわちベクトルのベクトル）で与えられているとしよう。未知のパラメータベクトル $b は PDL::Matrix オブジェクトでベクトルである。また共分散行列を multipiddle（ただし行列のベクトル）で与える。
問題を解くために本モジュールを使う前に、まず座って未知のパラメータに関する方程式を計算しなくてはならない。ガウス＝ヘルマートモデルでは、$y が観測値（ベクトルのベクトルである PDL::Matrix 上にあるものとする）で、$b が未知のパラメータベクトルであるような w($y,$b) = 0 という書式の方程式を指定する必要がある。この方程式 w($y,$b) = 0 はいわゆる「観測方程式（observation equation）」と呼ばれる。一般にはこれは「ベクトル方程式」すなわちスカラ方程式の集合であることに留意すること。
それから未知のパラメータベクトル $b も対して導入したい制限式について考える必要がある。（ユークリッド）距離を１にしたい、またはある値に対し１つのパラメータを設定したいとすれば h($b) = 0 という形式の方程式（ベクトル方程式）を指定する必要がある。
モデル w($y,$b)=0 及び h($b) = 0 に従いこれを最小化するためには観測関数 w($y,$b) 及び制限関数 h($b) のヤコビアンを指定する必要がある。観測ベクトル $y 及び未知のベクトル $b に従いヤコビアンが必要である。この段階では数値微分−ヤコビアンが不要−は不可能である。
モデルのコーディング
異なる n 個の観測ベクトル配列 @y があり、それぞれが PDL::Matrix ベクトルであったとする。たとえば PDL::Matrix のコンストラクタ vpdl で生成したものとする。推計にかけるためには $y = PDL::cat(@y) を用いて multipiddle $y を線上におく必要がある。
配列 @y の要素に対し共分散行列 @Sigma_yy の集合を適用する必要がある。
観測関数及びそのヤコビアン（未知のパラメータに関する制限式及びそのヤコビアンも）閉包として定義される。ヤコビアン閉包は PDL::Matrix 行列を返す。次元は @y 及び $b のベクトルの次元に依存する。
詳しくは t/ ディレクトリ及び example/ ディレクトリの例を見よ。
さらに推計のために以下のオプションを指定することができる。
verbose() ; verbose(2) ; verbose("/path/to/my/logfile")
変数が 1 に設定されていれば STDOUT に短い出力を行う。2 が指定されていれば詳細な出力がファイル "./GaussHelmert.log" に書き出される。これに替わるファイル名を指定することもできる。
verbose() の値のデフォルトは 0 、しなわち何も出力しない、である。
max_iteration() ; max_iteration(10)
条件が満たされない場合に計算を終了させる計算回数の最大値を取得または設定する。デフォルトは 15 である。
eps() ; eps($small_number)
繰り返し時の終了条件に使われる小さい値を取得または設定する。未知のベクトルをその標準偏差で割った値を最大誤差として参照される。デフォルト値は 10^(-2) であり、これは標準偏差の 1% の精度であることを意味する。
推定結果の出力
さあ、開始して結果を取得しよう。
start() ; start(verbose => 2)
推定を開始してオブジェクト自身を返す。
estimated_unknown()
推定した未知のパラメータベクトルを PDL::Matrix オブジェクトとして返す。
sigma0_squared()
推定した共分散因子を返す。
covariance_unknown()
推定した共分散行列を共分散因子なしに返す。
estimated_observations()
近似した観測値を multi-PDL オブジェクトとして返す。
estimated_unknowns_iterations()
未知のパラメータベクトルのリストを繰り返しごとに返す。
number_of_iterations()
繰り返しの回数を返す。
t/ 及び example/ フォルダを見よ。
観測値のベクトルのブロックを扱うサブクラス Statistics::GaussHelmertBlocks があるが。これはまだドキュメント化できていなくて、Math::UncertainGeometry を利用している。
今のところないが、冗長フラグで出力をチェックすることはできる。
おそらくあるだろう。しかし本関数は Math::UncertainGeometry ライブラリとともにかなりテストを通過している。
本モジュールは Perl Data Language を使用している。特に行列の演算のために PDL::Matrix を利用している。
２次元や３次元の予測点・線・平面を用いた本モジュールの使い方に関する複雑な例は Math::UncertainGeometry::Estimation を見よ。
Stephan Heuel (perl@heuel.org)
Copyright (C) 2000/2001 Stephan Heuel and Institute for Photogrammetry, University of Bonn. All rights reserved.
無保証である。本ソフトウェア／ドキュメントはいくかの条件のもとで再配布することができる。GNU Public License を見よ。GPL に関する詳細は http://www.gnu.org/copyleft/gpl.html
|
|
Updated : 2007/10/27 |