最小二乗法とロバスト推定

Size: px

Start display at page:

Download "最小二乗法とロバスト推定"

なおたておか
5 years ago
Views:

1 はじめに最小二乗法とロバスト推定 (M 推定 ) Maplesoft / サイバネットシステム ( 株 ) 最小二乗法はデータフィッティングをはじめとしてデータ解析ではもっともよく用いられる手法のひとつです Maple では CurveFitting パッケージの LeastSquares コマンドや Statistics パッケージの Fit コマンド NonlinearFit コマンドなどを用いてデータに適合する数式モデルを求めることが可能です一方で最小二乗法はデータ中にノイズをはじめとした例外値が含まれている場合に数式モデルが大きくズレてしまうことも多々ありますこの資料ではデータ中に含まれる誤差の影響を可能な限り少なくするための手法としてよく知られているロバスト推定 (M 推定 ) 法とその Maple 上での簡易的な利用法について解説しますまず最初に例外値が入ったデータでは最小二乗法によるフィッティングがどの程度ズレてくるのかを実際に確認してみましょう例外値の影響ここでは簡便に線形式 model = a$x Cb を考えますまた係数 a, b は次の通りとします : model d x/a$x Cb a d : b d 0.78 : model := x/a x Cb (.1) データのリストを変数 data に用意しますまず元となる X 座標 ( Xdata 変数で定義 ) Y 座標 ( Ydata 変数で定義 ) のリストを用意します Xdata d seq x, x = , 0.5 Xdata := 1.0, 1.5,.0,.5,.0,.5,.0,.5, 5.0, 5.5, 6.0, 6.5, 7.0, 7.5, 8.0, 8.5, 9.0, 9.5, 10.0 Ydata d seq model x, x = Xdata Ydata := 1.50, 1.805,.80,.8155,.00,.7905,.780,.7655, 5.50, 5.705, 6.80, , 7.00, , , , 9.150, 9.605, (.) (.) さらに Y 座標値には平均 0.01, 分散 0. の正規分布によるノイズを乗せますこれらのホワイトノイズは Statistics パッケージのコマンドを用いて生成します with Statistics : noise d convert Sample RandomVariable Normal 0.01, 0., nops Ydata, list noise := , , , , K , , , , K , (.)

2 , , K , , K , , K , , K , K Ydata d Ydata Cnoise Ydata := , , , , , , , , , , , , , , , , , , (.5) 作ったデータをグラフ描画して確認してみましょう with plots : gdata d pointplot Xdata, Ydata : display gdata このデータを元にした線形モデルへのフィッティングは Statistics パッケージの LinearFit コマンドで求まります fit1 d LinearFit 1, x, Xdata, Ydata, x fit1 := C x (.6) データ点のグラフと共に描画してみましょう gfit1 d plot fit1, x = 1..11, color = red :

3 display gdata, gfit x ここでもしもデータ点の Y 座標値にいくつかの大きな例外値が含まれた場合を考えてみましょう以下では Y 座標値の 5 番目と 18 番目にズレた値を設定してみます Ydata 5 d 7.8 Ydata 5 := 7.8 (.7) Ydata 18 d.08 Ydata 18 :=.08 (.8) もう一度グラフで確認してみます gnewdata := pointplot Xdata, Ydata : display gfit1, gnewdata

4 x ズレたデータに線形フィッティングを適用して数式モデルを確認してみます fit1 の結果と確認してみると傾きと切片の値も大きく変わったことがわかります fit d LinearFit 1, x, Xdata, Ydata, x fit := C x fit C x (.9) (.10) グラフでも確認してみましょう ( 元のフィッティングしたモデル式が赤色新しいデータに対するモデル式が緑色で描画されています ) gfit d plot fit, x = 1..11, color = green : display gnewdata, gfit1, gfit

5 x 次の章ではこのような例外値に対応するためのロバスト推定法を紹介します Tukey の Biweight 推定法前章で見たように最小二乗法はデータ点と数式モデルとの距離を最小化する手法です従って例外値のように大きくズレた値が含まれていた場合にはその影響を直接的に受けてしまうことを意味していますここで紹介する Tukey による Biweight 推定法は誤差の大きさに応じてデータ点に重みを付ける方法です具体的には誤差が大きい場合には重みを小さくすることで例外値の影響を小さくする手法ですそこで以下のような関数 biweight を定義します biweight d w, d /piecewise Kw % d and d % w, 1 K d w biweight := w, d /piecewise Kw % d and d % w,, 0 1 K d w, 0 (.1) ここで d : 誤差 w : 重み ( 誤差の許容範囲値 ) です例えば重み w =.0 とするとこの関数は次のようなグラフとなります : plot biweight.0, d, d = K5...5.

6 K K 0 d ここで定義した biweight 関数を元のデータとフィットしたモデルとの誤差値に適用して行きますまず fit 式による誤差値のデータを作りましょう errdata d Ydata K seq eval fit, x = xv, xv = Xdata errdata := K , K , K , K , , K , K , , K , , , , , , , , , K , (.) この誤差値に biweight 関数を適用して重み値のベクトルを作りますなおここでは重みの許容値を w =.0 とします weightvector d Vector seq biweight.0, di, di = errdata Vector column weightvector := Data Type: anything Storage: rectangular Order: Fortran_order (.) 計算された重み値のベクトルデータを用いて再び LinearFit コマンドに適用します重み値のベクトルは weights オプションに指定します

7 fit d LinearFit 1, x, Xdata, Ydata, x, 'weights'= weightvector fit := C x (.) 新しいモデル式 fit が得られましたグラフでその違いを確認してみましょう gfit d plot fit, x = 1..11, color = blue, legend = " 重みあり " : display gnewdata, gfit, gfit x 重みあり上記を見ると例外値の影響を少なくしたフィッティング結果が得られていることがわかりますまたフィッティングした数式モデルの結果を改めて確認すると次のようになっています 'fit1'= fit1; 'fit'= fit; 'fit'= fit; a, b fit1 = C x fit = C x fit = C x 0.975, 0.78 (.5) (.6) 次の章では Biweight 法を非線形のモデルフィッティングでも適用する方法を解説します

8 非線形モデルでの Biweight 法適用例線形の場合と同様の手順で非線形モデルフィッティングにも Biweight 法による M 推定を適用してみますまずデータの元となるモデル式を次のように定義します : restart model d t 1 p0 e p1 t p sin p t Kp Cp5 cos p6 t Kp7 Cp8 model := t/p0 e p1 t p sin p t Kp Cp5 cos p6 t Kp7 Cp8 (.1) 各パラメータ値は以下としますまたこのパラメータ値を代入したモデル式を f で定義します params d evalf p0 = 1.09, p1 = K0.89, p = 1., p = 5.08 p, p = K.1, p5 = 1.0, p6 = 1.87 p, p7 = K.8, p8 =. params := p0 = 1.09, p1 = K0.89, p = 1., p = , p = K.1, p5 = 1.0, p6 (.) = , p7 = K.8, p8 =. f d unapply eval model t, params, t f := t/1.09 e K0.89 t 1. sin t C.1 C1.0 cos t C.8 C. データ点の個数 N, 及びデータ点のリスト tvals, fvals を次のように用意します (.) N d 100 :.0 tvals d seq t, t = 0...0, N K1 fvals d seq f tv, tv = tvals : : 最後に正規分布のノイズ値を用意し fvals の値に加えます with Statistics : noise d convert Sample RandomVariable Normal 0.06, 0.09, N, list : fvals d fvalscnoise : グラフで確認してみます with plots : setoptions gridlines = true, axes = boxed : gmodel d plot f t, t = 0...0, color = red, legend = "Original Model" : gpoint d pointplot tvals, fvals : display gmodel, gpoint

9 1 0 1 t Original Model さてまずは Statistics パッケージの NonlinearFit コマンドでデータに対するモデルフィッティングを適用してみます fit1 d unapply NonlinearFit model t, tvals, fvals, t, t t fit1 := t/k e K sin t K C cos t K C (.) このフィッティングがどのような結果であるかをグラフで確認します gfit1 d plot fit1 t, t = 0...0, color = blue, legend = "NonlinearFit result" : display gmodel, gpoint, gfit1

10 1 0 1 t Original Model NonlinearFit result オリジナルのモデル式の振る舞いと比べると振幅や周波数で大きく異なっているのがわかりますそこで再び biweight 関数を定義し誤差値から重みデータを計算します誤差値を計算 : err d fvalsk seq fit1 t, t = tvals : biweight 関数を定義 : biweight d w, d /piecewise Kw % d and d % w, 1 K d w biweight := w, d /piecewise Kw % d and d % w,, 0 1 K d w, 0 (.5) 重みベクトルを計算 : weightvector d Vector seq biweight 9., d, d = err weightvector := Vector column Data Type: anything Storage: rectangular Order: Fortran_order (.6) 上で計算した重み値を与えて再度 NonlinearFit コマンドを適用してモデル式を求めてみます : fit d unapply NonlinearFit model t, tvals, fvals, t, 'weights'= weightvector, t (.7)

11 fit := t/ K e K t sin t C C cos t K C (.7) グラフで結果を確認してみます ( ここでフィッティングしたモデル式は緑色の曲線で示します ) gfit d plot fit t, t = 0...0, color = green, legend = "Biweight result" : display gmodel, gpoint, gfit1, gfit t Original Model NonlinearFit result Biweight result オリジナルのモデル式に極めて近い結果が得られました上記の例ではある特定の重み値を使いましたこの値が適切かどうかそしてより一般には重みの最適値はデータやモデル式によってまったく異なるため複数回の試行が必要になることに注意が必要です次の例ではこの w 値を探索するための簡易的なアプリケーションを用意していますロバスト推定を手軽に行うための使用例として参考にしてください重み値の対話的な変更によるフィッティング変化このアプリケーションを実行するにはまず下記のコード領域ボタンを押してコードを実行してくださいその後スライダーを動かすことで重み値を動的に変化させたフィッティング結果を確認することができます

weightvecfunc := proc(w) weight: 0 00 00 600 800 1000 ( 0.1) 9. 1 0 1 t Original Model Biweight result NonlinearFit result K0.05795869581 e K0.8880951009865 t 1.

関数の定義域を制限する関数のコマンドを入力バーに打つことにより関数の定義域を制限することが出来ます Function[ < 関数 >, <x の開始値 >, <x の終了値 > ] 例えば f(x) = x 2 2x + 1 ( 1 < x < 4) のグラフを描くには Function[ x^

関数の定義域を制限する関数のコマンドを入力バーに打つことにより関数の定義域を制限することが出来ます Function[ < 関数 >, <x の開始値 >, <x の終了値 > ] 例えば f(x) = x 2 2x + 1 ( 1 < x < 4) のグラフを描くには Function[ x^ この節では GeoGebra を用いて関数のグラフを描画する基本事項を扱います画面下部にある入力バーから式を入力し後から書式設定により色や名前を整えることが出来ますグラフィックスビューによる作図は後の章で扱います 1.1 グラフの挿入関数のグラフは関数 y = f(x) を満たす (x, y) を座標とする全ての点を描くことです入力バーを用いれば関数を直接入力することが出来その関数のグラフを作図することが出来ます