最尤推定の実践：最適化ルーチンの効果的利用

Geoff Morgan, nAG Ltd

1. はじめに

nAGライブラリには、統計モデルを適合させるためのルーチンが含まれます。一般的なケースでは、これらのルーチンを直接使用できますが、より専門的な例では、ユーザーはnAG E04章の最適化ルーチンを使用する必要があります。統計モデルは、一連のパラメータを用いてデータ値を観測する確率を与えます。パラメータの値は、モデルがデータに最もよく適合するように、データから推定する必要があります。例えば、観測値がポアソン分布に従うと仮定すると、値 $y$ を観測する確率は以下のようになります：

$P (y) = \frac{e^{- μ} μ^{y}}{y!}$

ここで、パラメータ $μ$ は分布の平均です。平均は、別の観測変数 $x$ を通じて、以下のモデルで関連付けられる場合があります：

$μ = e^{β_{0} + β_{1} x}$

$y$ と $x$ の複数の観測値が与えられた場合、パラメータ $β_{0}$ と $β_{1}$ の推定値を見つける必要があります。このようなモデルは、ポアソン誤差を持つ一般化線形モデル用のnAGルーチンG02GCFを使用して適合させることができます。G02GCFは、最尤推定法を用いてポアソン分布に従うデータに対して一般化線形モデルを適合させるルーチンで、様々なリンク関数（指数、恒等、対数、平方根、逆数）をサポートし、オフセット項や重み付けも扱うことができます。また、モデルの適合度を評価するための逸脱度や、パラメータ推定値の標準誤差、分散共分散行列なども提供します。

しかし、モデルが以下のようである場合：

$μ = \frac{1}{β_{0} + β_{1} x}$

ユーザーはnAGライブラリのE04章に含まれる最適化ルーチンを使用してモデルを適合させる必要があります。E04章は関数の最小化または最大化のためのルーチンを提供しています。この場合、最尤推定のために対数尤度関数を最大化することになります。E04ルーチンを使用することで、標準的なnAGルーチンでは直接サポートされていない複雑なモデルでも、最尤推定を行うことができます。

2. 最尤法の理論

統計モデルのパラメータの推定量を見つけるためによく使用される基準は、最尤法です。最尤推定では、観測されたデータが最も起こりやすい結果となるようにパラメータが選択されます。このような推定量を見つけるには、まず尤度関数を見つける必要があります。尤度は、確率変数の関数としてではなく、パラメータの関数として考えられるサンプルの確率分布です。独立した観測値のサンプルの場合、尤度は各観測値の確率分布の積です。例えば、ポアソン分布の場合、尤度は以下のようになります：

$L (μ) = \prod_{i = 1}^{n} P (y_{i}) = \prod_{i = 1}^{n} \frac{e^{- μ_{i}} μ_{i}^{y_{i}}}{y_{i}!}$

多くの場合、対数尤度はより単純な関数です。ポアソンの例では：

$\log L (μ) = \sum_{i = 1}^{n} (y_{i} \log μ_{i} - μ_{i} - \log (y_{i}!)) = \sum_{i = 1}^{n} (y_{i} \log μ_{i} - μ_{i}) + const$

ここで、constはパラメータに依存しない定数です。正規分布からの独立サンプルの場合、対数尤度は二乗和関数となり、最尤推定量は最小二乗推定量と同じになります。

3. 最適化ルーチンの利用

尤度関数または対数尤度関数を見つけたら、適切なnAG E04ルーチンを使用して最大化し、パラメータ推定値を得ることができます。正規分布の場合、使用すべき専門的な非線形最小二乗ルーチンがあります。

4. 分散共分散行列の推定

パラメータ推定値に加えて、ユーザーは一般的にパラメータ推定値の分散共分散行列も必要とします。これらはパラメータに関する仮説を検定したり、パラメータの信頼区間を設定したりするのに使用できます。また、パラメータ推定値間の相関を調べることも推奨されます。高い相関は、異なるパラメータ化を検討すべきであることを示すか、それが不可能な場合は、相関のあるパラメータに対してのみ共同で推論を行うべきであることを示します。内部最大値で見つかった最尤推定量の場合、大規模サンプルでは、分散共分散行列は対数尤度の2次導関数行列の期待値の負の逆行列で近似できることが示されています。しかし、多くの場合、導関数の期待値の計算が難しいため、代わりに観測値が使用されます。つまり、ヘッセ行列の負の逆行列が使用されます。

5. E04ルーチンの出力に関する注意点

多くのnAG E04ルーチンは、ヘッセ行列またはヘッセ行列のCholesky因子の近似を返し、多くのユーザーはこの情報を分散共分散計算に使用しようとして誤解しています。実際、返される情報は最適化プロセスを再開する文脈では問題ありませんが、一般的に統計計算には十分ではありません。特に、ルーチンがIFAILに非ゼロ値を返す場合（例えば、E04WDFのIFAIL=3）、ルーチンは通常、単位行列をヘッセ行列の近似として返します。

6. 推奨される手法

代わりに、ユーザーはヘッセ行列をより確実に推定する手段として、ライブラリに含まれているE04XAFを使用することをお勧めします。E04XAFは有限差分法を用いて勾配ベクトルとヘッセ行列の近似を計算します。このルーチンは自動的に適切なステップサイズを選択し、高精度な近似を提供します。勾配のみ、ヘッセ行列のみ、または両方を計算するオプションがあり、様々な最適化問題に柔軟に対応できます。E04XAFを使用することで、最尤推定のパラメータの分散共分散行列をより正確に計算できます。

また、ユーザーは、無制約非線形最小二乗ルーチンの1つを使用した後に分散共分散行列の推定値を提供するためにE04YCFを利用できます。E04YCFは、非線形最小二乗問題の推定回帰係数の分散共分散行列を計算するルーチンです。このルーチンは、e04fcf、e04fyf、e04gbfなどの非線形最小二乗ルーチンの後に使用されます。E04YCFは、解におけるヤコビアンを用いて分散共分散行列を推定し、パラメータの不確実性や相関を評価するのに役立ちます。ヤコビアンが特異な場合でも計算が可能です。このルーチンを使用することで、複雑な非線形モデルでも推定パラメータの統計的性質を適切に評価できます。