2012年8月31日金曜日

cos関数の線形最小二乗法

目的

データ(xi,yi)が複数(i=1n)与えられたとき、次のモデルのパラメータ(a,b,c)を線形最小二乗法により最適化する。

y=acos(x+b)+c

方法

  1. 加法定理とa=acosbb=asinbで変形し、3変数連立一次方程式の形にする。
  2. 平均自乗誤差を最小化する問題として線形最小二乗法の問題として解く。
  3. 求まった最適な(a,b,c)からa=a2+b2b=tan1(b/a)として元の変数を求める。

ソースコード

Jupyter NotebookのipynbファイルをGitHubで公開している.

STEP 1: 加法定理と変数変換でパラメータの線形和に

実際には、データ(x,y)には観測誤差が含まれて定式の等号は成り立たないので観測誤差をまとめてϵiと表すと、i番目のデータ(xi,yi)の関係は以下のように表現できる。

ϵi=acos(xi+b)+cyi

ここでの最適化とはこの全てのデータに対してϵiが小さくなる(a,b,c)を探すことなので、その指標として平均自乗誤差E=inϵi2/nを最小化する。加法定理で変形し、a=acosbb=asinbと置くと以下のように変形できる。

ϵi=acos(xi+b)+cyi=a{cos(xi)cos(b)sin(xi)sin(b)}+cyi=acos(xi)bsin(xi) +cyi

STEP 2: 誤差関数を微分する

今、平均自乗誤差Eを最小化しようとしているので、最適な値(a,b,c)ではEが極値になることが必要条件である。

Ea=0, Eb=0, Ec=0

以下、微分に関係ない定数nを左辺に移項して上式を実際に計算する。

nEa=ainϵi2=inaϵi2 =in2ϵiaϵi=in2{acos(xi)bsin(xi) +cyi}cos(xi)=2{aincos2(xi)binsin(xi)cos(xi)+cincos(xi)inyicos(xi)}=0nEb=in2{acos(xi)bsin(xi) +cyi}sin(xi)=0=2{aincos(xi)sin(xi)binsin2(xi)+cinsin(xi)inyisin(xi)}=0nEc=in2{acos(xi)bsin(xi) +cyi}=0=2{aincos(xi)binsin(xi)+cninyi}=0 

上記3つの方程式は、パラメータ(a,b,c)のみに注目すれば3変数の連立一次方程式である。大学生の線形代数の知識を使って、以下のように行列の式にすれば簡単に解ける。

[incos2(xi)insin(xi)cos(xi)incos(xi)incos(xi)sin(xi)insin2(xi)insin(xi)incos(xi)insin(xi)n][abc]=[inyicos(xi)inyisin(xi)inyi]

STEP 3: 変数変換を戻して答えを求める

求まった最適な(a,b,c)からa=a2+b2b=tan1(b/a)として元の変数を求める。プログラム上では関数atan2を使うとa=0の場合分けの必要がない。使い方は以下のとおり(ただし、a_=a、b_=b)。

b = atan2(a_, b_)