2002/10/26 梶本裕之

この文章は，最小二乗法，フィッティング，パラメータ推定の正しい解説テキストではありませんし，私もそれらの専門家ではありません．文中には私の誤解による部分が多々あるかと思いますので，くれぐれも信じないようお願いします．ただ私自身を含めて周りを見回してみると，制御や信号処理の割と高度なところまで勉強したはずなのに，いざ簡単なフィッティングになるとこっそりExcelに頼る他なかったり，インピーダンス同定で手も足も出ないことに気づく，という場面が多いようです．そういう状態の真っ只中にいる場合は参考になるかもしれません．

前半では最小二乗法と擬似逆行列を導入します．メインとなる後半では特にエンジニアリングでの応用例を幾つか挙げます．例といっても，各方面からバランスよく取ってきたものではなく身の回りに落ちていたものですので偏っています．

1. 最小二乗法

N個の入力()とN個の未知パラメータ()の積和によって1個の出力()が決まるシステム（線形システム）を考えます．

入力を変化させ， M回の測定を行ったとします．

結局

ただし

システムの未知パラメータxを求めるのが目的です．

もしM=Nなら行列Aは正方行列となり，逆行列が存在するので

でおわり．

となるとN回測定すればよいように思えますが，実際は測定には必ず測定誤差が存在しますから，誤差の影響を低減するために測定回数を増やします．このときM>Nとなり，Aは正方行列ではなくなり，逆行列は存在しなくなります．

1.1. 擬似逆行列

測定された出力ベクトルが，ある誤差()を含んだものであると仮定します．（Tは転置を表す）

もし自分の測定に自信があり，誤差には偏りが無いと思えるなら，誤差の「大きさ」が最小となるようなxをもっともらしいxとして受け入れることが出来るでしょう．この誤差の「大きさ」として二乗ノルムを採用するのが最小二乗法です．これは次のように計算されます．

*は複素共役転置を表し，実数であれば転置です．

この誤差の二乗ノルムの極値（がこの場合極小値であり，最小値でもある）を与えるxが，求めるパラメータベクトルです．

極値を与えるxを得るためにxに関して微分して

ここでをの擬似逆行列(Pseudo Inverse)と呼びます．今後と表記します．

もちろんもしシステムになんらかの情報，制約があると，より良い解が存在するはず，と予想できるでしょう．ここから先はとてもとても広い｢逆問題｣の分野が広がっているそうですが，よく知らないのでこれでおわり．

ただ話の肝は，が測定によって作られることと，その測定はあくまで自分で行うということ，良い（解が誤差の影響を受けにくい）を得るも得ないも測定次第，ということです．次から挙げる最小二乗法を用いたパラメータ推定の例は多くがセンシング関係ですが，センサは逆問題発生器であり，良いセンサは逆問題が解きやすいセンサであり，条件の厳しい逆問題を解くことに耽溺するよりも楽に解ける逆問題に直すための計測手段を探るべき，というのが正しい姿勢であることは間違いありません．

2. パラメータ推定の例

2.1. 力センサの較正（直線近似）

某社のフィルム状力センサは感圧インクを使用しており，表面に加えられた力(100gf程度)に応じて直線状にコンダクタンス（抵抗の逆数）が変化します．ただし比例ではなく，あるオフセットが存在します．コンダクタンス測定回路により電圧出力を得るシステムを作りました．

図 1 Vt=-5.00V,Rf=100kΩ,Rs：力センサ（注意：この回路のままだと力がかかっていない時Rs＝∞で熱雑音測定器になってしまうのでこの後にLPFを入れる必要あり．）

出力電圧V(volt)と表面応力F(gf)の関係は次のように書かれます．

ただしは比例係数，はオフセットを表すモデルパラメータです．を求める操作こそが較正(Calibration)であり，モデルパラメータの推定です．較正によって初めて「力センサ」として使うことが出来る（電圧から力を逆算できる）ことになります．

分銅を5種類(10gf,20gf,50gf,100gf,200gf)用意し，力センサに載せ，出力電圧を測定しました．分銅1種類につき10回の測定を行い，分銅を載せない場合を含め計60回のデータを得ました．

表 1重量と出力電圧

結局，

これは

の形ですから，擬似逆行列を用いて未知パラメータ，を求めることが出来ます．我々が行うべき操作は測定によってyとAを得ることであり，後は適当な行列ライブラリに任せればよいことになります．逆行列さえ計算できるライブラリを（そこら中に落ちています）持っているなら擬似逆行列は定義に従って求められます．Matlabであれば擬似逆行列はpinv(A)です．下の図が線形近似結果です（といいつつExcelを使ってますが…）

図 2センサ特性の線形近似結果

おまけとして，行列演算を要素ごとの計算で行い，を求めてみます．

結局，

最後の式は実験データの直線近似法として紹介される事が多いものです．特に学生実験レベルの化学実験や統計学の始めのほうでは直線近似ばかりやるので，最小二乗法イコール直線近似という誤解を生む原因となっているかもしれません．たしか初めてポケコン（関数電卓＋Basicみたいなもの）で作ったプログラムはこれでした．

また例では一種類の分銅に関して10回測定し，全部で50回のデータを取得，50組のデータ系列として処理していますが，あらかじめ一種類の分銅に関する10回の測定の平均を取ってから， 5種類の分銅による5回のデータとして計算しても結果は全く変わりません（説明省く）．一種類の分銅に関して多数回測定するのは誤差を小さくするという意味しかありません．

2.2. 多項式近似

データを直線近似してみるとどうも上手く合わないが，適当なモデルがわからない，という場面ではとりあえず多項式で近似してみることになります．多項式での近似はまず最大次数を決めるところから始まります．最大次数は直線の場合は１ですので，これを２，３と上げていき，データにフィットするかどうか見ていきます．ここでは２次の場合を見ます．

モデルは

から

に変わります．やはりデータをたくさん取ると

結局，

これは

の形ですから，擬似逆行列を用いて未知パラメータ，，を求めることが出来ます．実用的な状況では，モデルのパラメータを求めること自体は重要でなく，逆関数を求めること（ここでは電圧から力を逆算できるように較正すること）のみが要求されることも多々あります．そういった状況ではモデルの次数を適当に上げていくこの方法が有効になります．ただ，２次式までは逆関数も簡単ですが，３次以降どんどん難しくなることは覚悟すべきです．

2.3. フォトリフレクタを使った測距計の較正

図 3　フォトリフレクタ：三洋SPI-315の概観と内部概略（右の図はカタログより）

上の図は三洋製フォトリフレクタSPI-315です．赤外LEDとフォトトランジスタの組み合わされたものです．サイズは4mmx4mm位で1mm～10mmまでの反射板の距離を測ることが出来ます．全く同じサイズのものが各社から出ています．2mm～5mm程度の移動を高速に，でもアバウトに測る要求が出たので，距離と出力電圧の関係を測定により定式化する必要がありました．

図 4　測定回路と測定状況

上の図のような回路でLEDを光らせ，フォトトランジスタに生じる光電流を抵抗で電圧に変換して出力させました．回路のパラメータは個人的な制約が多々あった中での適当な値なので真似しないほうが賢明です．（このままだと電圧が高すぎるのでだんだん暖かくなり特性が変化していきます…出力が飽和しないように選べば大抵問題ないでしょう）

ハイトゲージの底部にフォトリフレクタを取り付け，反射板（ここでは鏡面反射するアルミ板）を取り付け，フォトリフレクタと反射板の距離を変化させながら出力電圧を測定していきます．

めんどくさくなってきたので昔のノートをそのまま載せます．

まずモデルを理論的に求めます．

次にデータのフィッティングを行います．フィッティングは最小二乗法で行いますが，非線形のものを無理やり線形にしていくため，最小二乗の本来の意義は失われていきます．こうした操作は最小二乗法をフィッティングに用いるためのテクニックとして必須なのですが，とんでもないフィッティング結果を出すことになりかねない諸刃の剣でもあります．

図 5　実測データ（○）と最小二乗フィッティング（実線）と逆関数によって電圧から距離を求めた場合の誤差．

上の図が実測データ（○）と最小二乗フィッティング（実線）と逆関数によって電圧から距離を求めた場合の誤差になります．誤差が50μm程度に抑えられていることが解ります．

このセンサを使う場合，グラフから明らかなように出力電圧に二値性があり，電圧だけからは距離が0.7mm以下なのか，以上なのか解りません．このためストッパーなどを設けて1mm～2mm以上

離して使うことになります．距離が離れると今回のモデルのように細かな話は必要なく，「出力電圧は距離の二乗に反比例」という簡易モデルで充分になります．今のモデルでは逆関数を求めるために３次方程式の解が必要でしたが，簡易モデルではずっと簡単になります．これを次に見てみます．

モデルは

ただしy：出力電圧[V]，x:距離[mm]，k1：比例項，k2：距離のオフセット．

すると

計測によりデータ群，を得ると，

となります，これは

の形であり，擬似逆行列を用いて未知パラメータk1,k2を求めることができます．

逆関数は

なので，この式により出力電圧yから距離xを求めることが出来ます．

下の図が新たな計測によるフィッティング結果です（距離はストッパーの位置を０としています）．逆関数によって電圧から距離を求めた場合の誤差も示しています．50μm程度の精度が出ていることが解ります．

図 6実測データ（○）と最小二乗フィッティング（実線）と逆関数によって電圧から距離を求めた場合の誤差（距離の二乗に反比例するとした場合）

2.4. 直交検波によるsin波の位相･振幅検出

後の章で述べるようなシステム同定を行う際，システムにsin波を入力し，出力のsin波の振幅，位相を測定することによって，伝達関数を求めることがあります．また例えば超音波やレーザの反射による距離計測では，反射波の時間遅れを精度よく計測するために反射波の位相遅れを測定します．このようにsin波の位相，振幅を計測する必要がある状況で多く使われる方法が直交検波です．まず前半で直交検波の原理を説明し，後半で直交検波と最小二乗法の関連について述べます．

直交検波の仕組み

問題：なるデータが与えられたとき，振幅Aと位相遅れθを求めたい（周波数ωは既知）．

まずsin(ωt)を掛け，1周期の時間区間Tで積分します．

積分区間が一周期であるため第２項が消えています．積分区間はN周期でも構いません．また適当な，1周期より充分長い時間であれば，第一項は累積していくのに第2項は累積しないため，相対的に第2項を無視することが出来ます．

同様にcos(ωt)を乗算し，1周期の時間区間Tで積分します．

二つの結果から振幅Aと位相遅れθは

として求めることが出来ます．ただしatan2はCやMatlab等大抵の数値計算ライブラリで定義されている関数で，の二値性を防ぎます．

直交検波と最小二乗法

直交検波は簡単なアナログ回路で実装することが出来ます（乗算回路と積分回路）．前述の超音波測距センサのように高速，実時間性が要求される場面ではアナログ回路で実装されます．しかしシステム同定のようにオフライン処理で充分な場合，また周波数ωを変化させる必要がある場合にはDAボードでsin波を入力し，システムの応答をADボードで取得，コンピュータ上で直交検波の演算を離散的に行うことになります．もちろんトリガによりDA，ADボードの同期をとることは大前提です．このとき離散的な直交検波が最小二乗法と等価になっている事を示します．

得られるデータは離散的な時系列データです．データ取得時刻を，データ系列をとします．

展開して

これは

の形であるため，擬似逆行列を用いて未知パラメータ，を求め，振幅A，位相遅れθを得ることが出来ます．

これを手作業で行ってみましょう．

Nを上手く選べば，は0に出来ます（時刻が等間隔の場合は一周期分だけ切り出せばOK）．この式がデータ系列に対する直交検波そのものであることは明らかでしょう(の意味に注目)．時刻は等間隔でなくてもかまいません．

2.5. インピーダンスの同定（電気抵抗とコンデンサ容量の同定）

上図のように抵抗とコンデンサによる回路があります．何らかの制約により，それぞれの抵抗値，コンデンサ容量を単独で測定することが出来ず， V(t)，I(t)の関係からR1,R2,Cを推定しなければなりません．このような場面は電気回路に限らず，音（音響インピーダンス），ロボット（機械インピーダンス）等，インピーダンスの関係するところに必ず出てきます．

I(t)を入力，V(t)を出力とします（逆でも構いません．両方測定できることが大切）． V(s)，I(s)をV(t)， I(t)のラプラス変換，Z(s)を複素インピーダンスとすると，

伝達関数は，

このが求まればR1,R2,Cが計算できる事は明らかでしょう．

測定は次の手順で行います．まず入力I(t)としてを加え、定常状態になったところで出力を計測します．Sin波入力時，定常状態では入出力比と位相差がそれぞれ伝達関数の大きさと偏角に対応し，

すなわち

となることが知られています（適当な古典制御の本参照，あるいはhttp://www.teu.ac.jp/kougi/ohyama/lecture/cii01/lec3.html等）から，振幅比，位相差を前述の直交検波等を用いて測定することにより，を得ることが出来できます．（定常状態であることが前提なので多少長めにsin波を入力し，最初を除いて切り出します）

ところで

ですが，この式を変形して，未知パラメータに関して線形な表現を求めておきます．

こうしてに関して線形な式を求めました．

Sin波の周波数を変化させ，それぞれの周波数における複素インピーダンスを求めます．測定周波数を，得られた複素インピーダンスのデータをとして，

これは

の形ですから，擬似逆行列を用いて未知パラメータを求め，R1,R2,Cを計算できます．

ここで述べた方法は周波数応答に関するフィッティングを行っており，視覚的には実測したボード線図に上手く合うモデルを求めるものです．幾つかの周波数のsin波を入れなければならない点で多少時間がかかります．あと，教科書的にはインパルス入力に対する出力をラプラス変換してやればそのまま伝達関数になる，とありますが，どうやってインパルスを作って入れるの？という所で固まってしまいます．というか，そんな入力を入れられたら壊れるのでは？より重要なことは，システムを線形と見なせる領域は必ずしも広くないので，｢普段使う周波数の範囲｣かつ｢普段使う振幅｣の入力に対する応答を測るべきということです．インパルスという特殊なものは相当そのシステムの線形性を信じていなければ使えないもので，結局ここで紹介した古典的な方法は割と有効ということになります．

時間軸で議論することも出来ます．以下はロボットアームの動力学パラメータ（ヤコビアン）同定等によく使われる方法です．運動方程式の形のままで同定します．

sが微分演算子であることを思い出すと

結局，時刻での，，，さえ求まれば

これは

の形ですから，擬似逆行列を用いて未知パラメータを求め，R1,R2,Cを計算できます．

もちろんV(t),I(t)の微分はV(t),I(t)の差分をとって求めれば良いです．

この方法の良いところは，特別な入力を与えなくても実際に使うための入力で同時に同定もできるところです．よく言われる落とし穴は微分を取るところでノイズの影響が大きくなるということです．一般には微分項も直接計測できればそれに越したことは無いとされています．例えばモータならポテンショメータだけでなくタコメータを使って微分項を直接取ってやるとか．または各項を時間積分した形に直して微分が現われないようにするとか，各項のノイズの大きさを評価しておいて重みつきの最小二乗法にするとか，普段使う周波数範囲がまんべんなく入るようにチャープ波を入れるとか，きりがありません．ただ我々は同定問題の専門家では無く，そのシステムを同定する目的が何であるか知っているわけですから，その目的が達成できる程度に使えるならどんな方法でも良く，むしろ大上段な同定問題の教科書を使うよりは最も単純な方法から試していくべき，と個人的には思います．

ここから先はこの文章の守備範囲ではありません．キーワードは多分状態方程式，ARモデル，カルマンフィルタ，etc…？この先はいくらでも本に載っているようです．