ロバスト推定法(Tukey’s biweight)

ロバストとは?

ロバスト(Robust)という言葉を辞書で調べると、「頑健なさま、がっしりした様子」という意味が載っています。

 

画像処理的には、ノイズや影、明るさの変動などの影響を受けにくく、安定した処理結果が得られる、ぐらいの意味として捉えれば良いでしょう。

 

ロバスト推定法

通常の最小二乗法では下図のように、数点の大きな誤差が含まれるだけでも、近似した直線が大きくズレてしまう場合があります。
この誤差の影響をできるだけ受けないようにしたのが、ロバスト推定法です。

 

処理アルゴリズム

まずは最初に、通常の最小二乗法を行います。
(今回は簡単にするため一次式で近似する例で紹介します。)

 

直線の式(y=ax+b)で近似する場合は以下のように行列を用いて解くことができます。

 

この近似した直線から遠く離れたデータを除去するだけでも、大きな誤差のデータの影響を受けなくなりそうですが、ここでは、誤差の大きさに応じて重みを付けるTukeyのBiweight推定法という手法を紹介します。

 

近似データ(Xi、Yi)と近似直線との誤差 d=Yi – (aXi + b) を用いて、誤差が大きければ大きいほど、
最小二乗に与える影響力(重み)が小さくなるように、以下のような式を用いて重み計算します。

 

d < -Wの場合

-W <= d <= Wの場合

W < dの場合

 

ただし、大文字のは誤差の許容範囲を示します。

 

この重みの関数w(d)をグラフで示すと、このようになります。

 

 

上図を見ても分かるように誤差の絶対値がWを超えると重みが0となり、次に行う重み付き最小二乗法では近似に影響を与えなくなり、逆に誤差が0に近づく程、重みが増し、近似への影響が大きくなります。

 

この重み wi を各近似データ(Xi、Yi)に関して計算し、wiを付加した最小二乗法を再度行います。

 

 

こうして求めた近似式(y = a’X + b’)は最初に求めた近似式よりも、測定データに近づきます。
この処理を誤差の許容範囲()を小さくしながら、誤差が少なくなるまで繰り返すことで、より測定データに近づいた近似式を得ることができます。

 

今回の例では一次式について示していますが、他のn次式などについても、同様に行列の各要素に重みを付加することで、ロバスト推定を行うことができます。

 

参考文献

 

使える数学へ戻る

 

最小二乗法の最適化(高速化)

最小二乗法を用いて曲線近似するとき、近似するX座標の値が等間隔に並んでいる場合、奇関数、偶関数の特性を使って最小二乗法を高速に処理することが可能になります。

 

例えば、画像中の線の中心位置や、エッジの位置をサブピクセル単位で求める場合などに有効です。

 

 

以下からは線の輝度値を2次曲線で近似して、線の中心を求める例を示します。

 

X方向の輝度値が最大となる座標付近の輝度値を抜き出し、2次式で近似すると
X座標が180と181との間に最大値がありそうな事が分かります。

 

 

この輝度値を二次式(Y = aX2 + bX + c)で近似する場合、下記の行列の式を解くと二次式が求まります。

 

 

ここで、逆行列の中は奇関数、偶関数の集まりなので、輝度値データのX座標を中心座標(X = 180.5)で割り振ると

 

 

となります。
よって、奇関数、偶関数の特性から、行列の式は

 

 

となり、行列の計算量を減らすことができます。

 

また、近似するデータの個数が固定できるなら、逆行列の部分はあらかじめ計算しておく事ができるので、さらに計算量を減らすことができます。

 

これで、求める線の中心座標は近似曲線の頂点位置、つまり微分の値が0となる座標(Y’ = 0の座標)なので、

線の中心座標 = Xの中心座標 – b / 2a
(今回の例ではXの中心座標 = 180.5)

 

となり、より高速に線の中心座標を求めることができるようになります。

 

この前提条件とした『X座標の値が等間隔に並んでいる場合』という部分ですが、画像処理においては、画素のX座標、Y座標を用いて、輝度値に関して近似する場合が多いので、for文などで、ある関数の総和などを求める場合には、この奇関数・偶関数が使えないか?検討してみると良いと思います。

 

このテクニックは計算量が格段に少なくなるので、効果的な場合が多いと思います。

 

使える数学へ戻る

 

最小二乗法をExcelで解く

エクセルではグラフの機能で近似曲線の追加というのはありますが、任意の式で近似したい場合や、近似式の係数を取得したい場合にはエクセルの近似曲線の機能だけでは物足りない場合があります。
そんなときにエクセルでまじめに最小二乗法で解く方法を紹介します。

 

まずは、最小二乗法を解くとき、逆行列を用いるところまではご理解願います。
詳しくはこちらで紹介しています。

 

最小二乗法
一般式による最小二乗法

 

また、Excelで行列を計算する方法はこちらを参照下さい。

 

行列の積、逆行列、転置行列の計算

 

以下では円の最小二乗法を例にとって説明をしたいと思います。

 

円の最小二乗法を解くときの円の一般式を

 

X2 + Y2 + AX + BY + C = 0

 

としたときの行列の式は

 

 

です。
まず、行列の式に使われているの値を求めます。
はエクセル関数の「SUM」関数を用いればよいので、簡単に求まります。

 

 

まず、行列の式の各値を各セルに書いておきます。
次に逆行列を求めるのですが、エクセル関数に「MINVERSE」というのがあります。

 

逆行列を表示したい複数のセル(下記の例では3×3のセル)を選択します。
次に

 

=MINVERSE(逆行列を求める複数のセル)

 

と入力します。ここでは文字を入力するだけでEnterキーは押さないこと

 

 

 

で、次が最大のポイント!
通常ではEnterキーを押したいところですが、行列の計算では

 

Ctrlキー  + Shiftキー + Enterキー

 

を同時に押します。
これで、逆行列の値が求まります。

 

 

あとは行列の積が計算できれば、求めたいA,B,Cの値を求めることができるのですが、
行列の積の計算にはエクセル関数の「MMULT」という関数を用います。

 

逆行列を計算したときと同様に行列の積の結果を表示したいセル複数のセルを選択し、

 

=MMULT(行列1,行列2)

 

を入力します。ここでも最後はEnterキーだけではなく

 

Ctrlキー  + Shiftキー + Enterキー

 

を押すので注意して下さい。

 

 

これで求めたかった未知数A,B,Cが求まります。

 

 

今回は円の近似なので、A,B,Cの値を用いで円の中心座標(a,b)、円の半径rを求めるには

 

a = -A / 2
b = -B / 2
r = SQRT(a2 + b2 – C)

 

なので、

 

 

となり、円の近似式を求めることができます。

 

このエクセルによる円の最小二乗法のファイルはこちら↓

LSM_Circle.xls

 

また、一般的なn次式による最小二乗法を解いたのファイルはこちら↓

n_LeastSquare.xls

 

使える数学へ戻る

 

一般式による最小二乗法(円の最小二乗法)

n次曲線による最小二乗法につていは説明しましたが、円や楕円、その他の一般式についても最小二乗法による近似は可能です。
今回は円の最小二乗法を例にとって説明しますが、他の一般式についても要領は同じです。

 

まず初めに近似する一般式を作成します。

円の場合、円の中心座標(a,b)、円の半径 r とすると

 

(X – a)2 + (Y – b)2 = r2 ・・・①

 

となります。

n次曲線による最小二乗法ではX座標がXiの時のY座標Yiと近似曲線上のY座標f(Xi)との差の二乗が最小になるように計算していましたが、円の場合、X座標がXiの時にY座標が存在しなかったり、Y座標が2点存在してしまう場合もあります。

 

そこで、今回はY座標の差の2乗の総和を求めるのではなく、①式を=0の式に変形し、
点の座標(Xi,Yi)を代入し、その式の2乗の総和を求めます。
(n次式近似でやった計算も結局は近似式=0になるように変形していた訳ですが...)

 

①式を=0になるように変形し、2乗の総和は

 

∑{(Xi – a)2 + (Yi – b)2 – r22 = 0 ・・・②

 

となります。

ただし、このまま未知数a,b,rに関して②式を偏微分してもa,b,rは4次関数となるので、各未知数の偏微分=0の時が最小とは限りません。

 

そこで、②式を展開し、下の式のように置きます。

 

∑{Xi2 + Yi2 + AXi + BYi + C}2 = 0 ・・・③

 

ただし

A = -2a ・・・④
B = -2b ・・・⑤
C =  a2 + b2 – r2 ・・・⑥

 

③式に関してA,B,Cについて偏微分すると

 

∂/∂A = A∑Xi2 + B∑XiYi + C∑Xi + ∑Xi3 + ∑XiYi2 = 0 ・・・⑦
∂/∂B = A∑XiYi + B∑Yi2 + C∑Yi + ∑Xi2Yi + ∑Yi3 = 0 ・・・⑧
∂/∂C = A∑Xi + B∑Yi + C∑1  + ∑Xi2 + ∑Yi2 = 0 ・・・⑨

 

⑦、⑧、⑨式を行列を用いて解くと

 

 

となり、逆行列を求めるとA,B,Cが求まります。

さらに④、⑤、⑥式より円の中心座標(a,b)、円の半径 rが求まります。

  \(a=-\frac{A}{2}\)

  \(b=-\frac{B}{2}\)

  \(r=\sqrt{a^{2}+b^{2}-C}\)

 

使える数学へ戻る

 

関連記事

最小二乗法の計算方法

最小二乗法の計算方法

最小自乗法という人もいますが、私は最小二乗法派です。

最小二乗法とは?

点の集まりから、近似直線(曲線)などの近似式を解くための手法です。
単に直線や二次曲線、三次曲線などに限らず、最小二乗法の応用範囲は広くあります。

最小二乗法は、言葉の通りに2乗を最小にする方法なのですが、この2乗と言うのが、方程式上の理論値と近似するデータとの差の2乗の合計を最小にする方法となります。と言っても分かりづらいので、具体的な最小二乗法の例を示します。

 

最小二乗法の計算方法

今回は説明を簡単にするため、f(X)= aX + b の直線でデータを近似する場合で説明します。

 

 

X = Xi のとき、直線上のY座標は

$$Y_{i}=f(X_{i})=aX_{i}+b$$

より、点と近似直線との差は

$$Y_{i}-f(X_{i})=Y_{i}-(aX_{i}+b)$$

となります。

 

この差の合計が一番小さい時にもっとも点の集まりを直線で近似できていることとなりますが、
そのままこの差の総和を計算すると差の値がプラスの時とマイナスの時があるため、例えば下図

 

 

のように、正の差と負の差が同じようにある場合、差の合計が小さくなってしまうため、差を2乗してから合計した値が最小になるようにします。

 

よって、差の2乗の合計は

$$\sum(Y_{i}-(aX_{i}+b))^{2}\\
=\sum(a^{2}X_i^2+2aX_{i}b+b^{2}-2aX_{i}Y_{i}-2bY_{i}+Y_i^2 )$$

・・・式①

となり、この総和の値が最小となるaとbを求めれば近似直線を求められることとなります。

 

この式を良くみると、aとbの二次関数でできているため、この二次関数が最小となる値を求めればよい。
二次関数が最小となるのは微分した値 = 0 の場合であるから、式①をaとbで偏微分した値が
0(ゼロ)となる値を求めればよい。

 

式①をaとbでそれぞれ偏微分すると

$$\frac{\partial}{\partial a}=2a\sum X_i^2+2b\sum X_{i}-2\sum X_{i}Y_{i}=0\\
\frac{\partial}{\partial b}=2a\sum X_i+2b\sum 1-2\sum Y_{i}=0$$

より

$$a\sum X_i^2+b\sum X_{i}-\sum X_{i}Y_{i}=0\\
a\sum X_i+b\sum 1-\sum Y_{i}=0$$

となり、未知数がaとbで式が2本あることから、aとbの連立法的式を解くことで未知数a,bを求めることができます。

 

連立方程式を行列で表すと

$$\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\left(\begin{array}{c}\sum X_i^2 & \sum X_{i}\\ \sum X_{i} & \sum 1\end{array}\right)^{-1}\left(\begin{array}{c}\sum X_{i}Y_{i}\\ \sum Y_{i}\end{array}\right)
\\
\left(\begin{array}{c}a\\ b\end{array}\right)=\frac{1}{\sum X_i^2\sum 1-\sum X_{i}\sum X_{i}}\left(\begin{array}{c}\sum 1 & -\sum x_i\\ -\sum x_i & \sum X_i^2\end{array}\right)\left(\begin{array}{c}\sum X_iY_i\\ \sum Y_i\end{array}\right)$$

となり、行列を解くと未知数a,bが求まり、近似直線を求めることができます。

 

ちなみに、二次式(aX2 + bX + c)で近似すると

$$\left(\begin{array}{c}a\\ b \\ c\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

三次式(aX3 + bX2 + cX + d)で近似すると

$$\left(\begin{array}{c}
a\\
b\\
c\\
d
\end{array}\right)=
\left(\begin{array}{c}
\sum X_i^6 & \sum X_i^5 & \sum X_i^4 & \sum X_i^3\\
\sum X_i^5 & \sum X_i^4 & \sum X_i^3 & \sum X_i^2\\
\sum X_i^4 & \sum X_i^3 & \sum X_i^2 & \sum X_i\\
\sum X_i^3 & \sum X_i^2 & \sum X_i & \sum 1
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum X_i^3Y_i\\
\sum X_i^2Y_i\\
\sum X_iY_i \\
\sum Y_i\end{array}\right)$$

と規則的に変化するため、n次式近似のプログラムも簡単に作ることができます。
ただし、行列の部分を二次元配列と考えると、近似する次数(n)が上がって行くと、配列の先頭に値が追加されていくのは、扱いにくいので、n次式近似のプログラムを作る場合は、近似式を

$$Y=a+bX+cX^2+dX^3 \cdot\cdot\cdot\cdot$$

として、行列の表現は

$$\left(\begin{array}{c}a\\ b\\ c\\d\\:
\end{array}\right)=
\left(\begin{array}{c}
\sum 1 & \sum X_i & \sum X_i^2 & \sum X_i^3 & .. \\
\sum X_i & \sum X_i^2 & \sum X_i^3 & \sum X_i^4 & .. \\
\sum X_i^2 & \sum X_i^3 & \sum X_i^4 & \sum X_i^5 & ..\\
\sum X_i^3 & \sum X_i^4 & \sum X_i^5 & \sum X_i^6 & ..\\
: & : & : & : & :
\end{array}\right)^{-1}
\left(\begin{array}{c}
\sum Y_i\\
\sum X_iY_i \\
\sum X_i^2Y_i\\
\sum X_i^3Y_i\\
:
\end{array}\right)$$

として、n次式で近似した最小二乗法を求めます。

 

この計算をエクセルで行ったものはExcelによる最小二乗近似からダウンロードできます。
逆行列を求めるプログラムにはガウスの消去法などが用いられます。

 

最小二乗法での注意点

三次式(a+bX+cX2+dX3)の近似を例にとってみると、近似処理をする過程でXに関する計算式はX~X6まで計算されます。

 

ここで、例えばXの値の範囲が

 

0.01~1000

 

だったとすると、Xに関する計算では

 

0.000000000001~1000000000000000000

 

までの計算がされます。
逆行列の計算ではこれらの値に関して足し算や掛け算がなされる訳ですからコンピュータの
計算誤差(桁落ち、情報落ちなど?)が発生する場合があります。
これを防ぐために、最小二乗法のプログラムでは近似式をそのまま使うのではなく、Xの値から
Xの平均値(X)を引いた計算式

$$Y=a+b(X-X_b)+c(X-X_b)^2+d(X-X_b)^3$$

で処理を行うのが一般的です。

 

最小二乗法の関連記事

一般式による最小二乗法(円の最小二乗法)

最小二乗法をExcelで解く

最小二乗法の最適化(高速化)

ロバスト推定法(Tukey’s biweight)

 

使える数学へ戻る