プロローグ
この記事では、SfMなどで使われるバンドル調整の根幹をなす、最小二乗法・最尤法の性質を扱います。
2変数X, Yに関するn個の観測データがあるとき、
xを説明変数とし、Yを目的変数とする回帰モデル：
Y = f(X) + ε
を考えて、関数fに含まれるパラメータを最小二乗法で推定する状況を考えます。
もし、Yの観測誤差を表すεが互いに独立で、平均0、標準偏差σの同一の正規分布に従うならば、最小二乗法は最尤法に一致します。すなわち、最小二乗法によって推定されたパラメータは、最尤推定量となります。
もしさらに、関数f(X)が線形つまりf(X) = a X + bの形ならば、最尤推定量は不偏性：
E(aの推定誤差) = E(bの推定誤差) = 0
を獲得します。つまり、推定のバイアス（平均的な誤差）がありません。
しかし関数f(X)が非線形の場合、一般には最尤推定量は不偏推定量ではなく、漸近不偏性を持つにすぎず、従ってnが十分に大きくなければバイアスが無視できません。このことを、Rを用いた数値実験で確認してみました。
数値実験の手順
① 回帰モデルを設定します。今回は次式を用い、aの真値は1.23, bの真値は0.456としました。
Y = a ／ (X^2 + b) + ε
② この回帰モデルと、εに関する上記の下線部の仮定に基づいて人工的なデータを生成します。今回は、1試行につき(X, Y)を30組作りました。具体的には、Xを区間[0, 1]での一様分布に従う互いに独立な疑似乱数で与え、対応するYの観測値は X, Yの真の関係式：
Y = 1.23 ／ (X^2 + 0.456) + ε
で計算しました。 ただし εは、平均0, 標準偏差0.2の正規分布に従う互いに独立な擬似乱数で与えました。
③ 最小二乗法でa, bを推定します。つまり、Yに関する30個の残差：
a ／ (X^2 + b)－Yの観測値
の二乗和（ないしRMS）が最小になるように、aとbを最適化（調整）します。今回は、準ニュートン法を応用した方法で最適化を行いました。
④ 上記②③を多数回繰り返し、aやbの推定誤差（推定値と真値との差）の統計をとります。今回は3000試行繰り返しました。
数値実験の結果：ある試行の例
まずイメージのため、ある試行の結果の例を示します。
ある試行におけるa, bの推定値は
a: 1.1766221, b: 0.4224552
、Yに関する30個の残差のRMSは
0.1582793
となりました。
また、Yの推定値（あてはめ値）と観測値の散布図は次のようになりました。赤線は直線 Y=Xです。
数値実験の結果：3000試行の統計
さて、3000試行（3000個の標本）に関する、a, bの推定値の標本平均は、
a: 1.2377946, b: 0.4600379
となりました。真値：
a: 1.2300000, b: 0.4560000
と比べると、a, bとも平均的には大きめに推定されていることになります。
a, bの推定「誤差」の標本平均に直すと、次のようになります。
a: 0.0077946, b: 0.0040379
これは統計的に有意なバイアスなのでしょうか？
3000試行に関する、a, bの推定値の不偏標準偏差は、
a: 0.1039258, b: 0.0471030
です。これを√3000で割って標本平均の標準誤差を求めると、
a: 0.0018974, b: 0.0008600
青字どうし、緑字どうしを比べると、推定誤差の標本平均は、その標準誤差の4倍を超えています。
正規分布に従う確率変数において、標本平均の絶対値が標準偏差の4倍を超えることは極めて稀（確率が0.006%未満）ですから、今回の実験で得られた推定誤差の標本平均は0と統計的に有意に異なり、最尤推定量に統計的に有意なバイアスがあると言えます。
おわりに
Yを目的変数とする回帰モデル：
Y = f(X) + ε
を例に、f(X)が非線形の場合、最尤推定量が不偏推定量ではないことを、1ケースだけですが数値実験で確認できました。
バンドル調整のモデル（連立された共線条件式からなるモデル）は回帰モデルではありませんが、最尤推定量が一般には不偏推定量でないということは重大な意味を持ちます。
画像に写ったタイポイントの観測誤差が「互いに独立で同じ正規分布に従う」という単純な条件下でも、推定されるカメラパラメータにバイアスが生じうるということだからです。
ただし、バンドル調整に関する文献によっては最尤推定量が不偏であると書いてあるので、バンドル調整に関しては例外なのかもしれません。このあたりは勉強して決着をつけたいところです（1年くらい放置状態）。
