« 電子を大きさのない点であると考える背景 | トップページ | 安倍政権に思う »

2006年9月22日 (金)

非線形最小二乗法(JEA式の作成過程)

 今日は窒素酸化物濃度などを計算する解析的な低煙源拡散式のパラメータの推定のために用いたことがある非線形最小二乗法について解説してみたいと思います。(直角風のJEA式)

 地上の道路上を走行する自動車からの窒素酸化物(NOx)などの汚染物質排出のソースを,有効煙源高さ(effective height):he=0 (m)のy軸上で単位長さ当たりの定常的な排出強度がQ(Nm3/(ms)の無限線煙源としてモデル化します。

 

 そして,道路に直角のx方向に風が吹いている設定で,風速も拡散係数も高さz(m)のベキ乗則に従って鉛直上方に向かって増加するというモデルを想定します。

 

 このときの拡散方程式の解である汚染物質濃度を,垂直距離x(m)と高さz(m)の関数と考えてC(x,z)で表わすと,これはAをQに比例するあるパラメータとして,C(x,z)=Ax-sexp(-Bzp/x)なる形の式(Robertsの式)で表わされる,ことがわかります。

 様々な環境の下で数回に渡って行なわれた実際の実験時の煙源長さはもちろん有限なのですが,とりあえず,これを無限大長さで近似します。

   

 

 実はガウスの誤差関数でもって有限長さの効果を取り入れることもできるのですが,今日の話題では割愛します。

 

 実験はトレーサーガスを地上にある有限長さの直線状のパイプに均等間隔で開けた穴から拡散させるというものです。

そのパイプ状の発生源をy軸としたとき,n個の観測点座標:(xi,zi)(i=1,2,..n)において,風向が直角に近い環境のときの実測濃度:ciと先に設定した低煙源拡散式による計算値:Ci=C(xi,zi)とを比較することによって,逆にその式のパラメータA,s,B,pを推定します。

 具体的には,C(x,z)=Ax-sexp(-Bzp/x)の右辺をf(A,s,B,p;x,z)と書いて,i=f(A,s,B,p;xi,zi)として,この計算値と実測値ciの誤差の二乗和:S(A,s,B,p)≡∑(Ci-ci)2を最小にするパラメータA,s,B,pを求めるわけです。

 

 これは,次に述べる非線形最小二乗法のパラメータ数が4個のときの例になっています。

 ここでは,より一般的に未知パラメータがr個あるとし,それらを順にA1, A2,..,Arとします。

 

 上のケースではr=4で,A1=A,A2=s,A3=B,A4=pです。

 

 そして,Ci=C(xi,zi)=f(A1,A2,..,Ar;xi,zi)とし,S(A1,A2,..,Ar)≡∑(Ci-ci)2と定義するわけです。

 誤差の二乗和Sが最小となる条件は,

 

 通常の線形回帰の計算の場合と同様,必要条件として

 i(A1,A2,..,Ar)≡∂S/∂Ai=2∑(Ck-ck)∂Ck/∂Ai=0 (i=1,2,..,r) で与えられます。

 

 そして,この r 個の連立方程式を解くことにより未知数A1,A2,..,Arを求めることが主目的となります。

 

 これらの方程式は一般に非線形ですから,こうした非線形回帰によるパラメータ推定の方法を非線形最小二乗法と呼ぶことにします。

 具体的には,i(A1,A2,..Ar)=0 (i=1,2,..r)を線形近似することにより多変数ニュートン法(Newtonian method)を実行します。

 

 そこで,まず連立方程式を線形近似します。

 

 すなわち,初期値としてAj=Aj0 適当に与えた後に,

 

 0 = i(A1,A2,..Ar)

  =i(A10,A20,..,Ar0)+∑(∂Fi/∂Aj)|A=A0(j-Aj0)

 

 と近似します。

これを行列形式で書くため,Dという行列をD=(dij)≡(∂Fi/∂Aj)で定義し,特にD0(dij0)≡(∂Fi/∂Aj)|A=A0とします。

 

一方,列ベクトルA≡t(A1,A2,..,Ar)で定義し,特に

0t(A10,A20,..,Ar0)とします。

 

さらに,Fi0i(A10,A20,..,Ar0)を成分とする同様な列ベクトルを0と定義します。0t(F10,F20,..,Fr0)です。

 

こうすれば,先の線形近似は 00+D0(0)と表わされます。

 

これを単純に解くと,00-10と近似されることになります。

ここで0-10の逆行列です。

得られた00-10を,改めて0として代入してD0-10を計算し,00-10収束するまで繰り返します。

 

すなわち,m+1mm-1mの漸化式において誤差|m+1m|の相対値が十分小さくなるまで繰り返し計算します。

 

具体的には,100-10,211-11,..,m+1mm-1mを用います。ただし,Dk(dijk)≡(∂Fi/∂Aj)|A=Akです。

 

m→ ∞では,|m+1m|→ 0,それ故m0 となることが予想されるため,この反復計算を実行するわけです。

 

こうすれば,m→ ∞ でのmt(A1m,A2m,..,Arm)の極限値としてi(A1,A2,..,Ar)=0 (i=1,2,...r)を満たすt(A1,A2,..,Ar)が得られるはずです。

さらに,実際には収束を速くするために加速係数:wを与えて,m+1m-wm-1mとする方法などを用いた方がいいと思います。

  

(※ここではwを加速係数であると称してはいますが,それは広い意味でありw>1を意味するわけではありません。実際,計算ではwとして2.0や逆に0.01など収束が悪い場合には,さまざまな値で試行錯誤しました。)

実際のr=4 の計算では,かつてFortranを使ってプログラミングしましたが,初期値を適切に取ると各環境のケースについて大体十数回の繰り返し計算でパラメータの推定値が得られました。

 

そして濃度の実測値とその推定パラメータによる計算値との相関係数としては,0.9 前後の値が得られ,回帰係数の値も0.8 から 1.2程度となって,仮定した計算式が良い近似になることがわかりました。

http://fphys.nifty.com/(ニフティ「物理フォーラム」サブマネージャー)                                       TOSHI 

http://blog.with2.net/link.php?269343(ブログ・ランキングの投票)↑ここをクリックすると投票したことになります。

|

« 電子を大きさのない点であると考える背景 | トップページ | 安倍政権に思う »

311 .数値計算・調和解析・離散数学」カテゴリの記事

コメント

 こんばんは。。hirotaさん。TOSHIです。コメントありがとうございます。

 なるほど線型な場合に導かれるニ乗和が最小になる条件を、非線型な場合にも適用すれば確かに2階微分は不要ですね。

 まあ、約30年くらい前に仕事上の必要にせまられて「アルゴリズム」を作る必要があったころは私も20代で社会人になったばかりで、そうした素養もなく「非線型最小二乗法」という言葉さえ知らなかったわけです。
 
 単にニ乗和が最小になる条件を求め、それをニュートン法で解くにはどうしたらいいかと考えた結果として「アルゴリズム」ができたわけで、何年か後にこれが「非線型最小二乗法」という方法と同じだとわかったわけです。

 そして当時、Acos6というNECと東芝の開発した大型コンピュータで1つの観測当たり30ポイントくらいのデータの数百観測についてFortranで4個か5個のパラメータを持つ、とても特異な形をした試行関数をいくつか設定して約2年間ほど計算をしましたが、形が特異であるにも関わらず、この2階微分する方法は結構うまく収束しました。

 変形ベッセル関数や楕円関数などは微分するだけでも関数近似式をつくらねばならないので苦労したのですが、コンピュータだとある関数の微分を計算するサブルーチンがあれば合成関数の微分は代入するだけなので2階程度の微分は机で計算するよりも比較的楽でしたね。

             TOSHI 

投稿: TOSHI | 2007年7月18日 (水) 19時45分

非線形最小二乗法の普通の計算法は、C_k を線形近似するだけです。
[C_k] を C_1 ~ C_n を要素とする列ベクトル、[A_i] を A_1 ~ A_r を要素とする列ベクトル、[∂C_k,i] を ∂Ck/∂Ai を要素とする行列として
  [C_k] = [C_k(0)] + [∂C_k,i] ( [A_i] - [A_i(0)] )
と線形近似し、| [C_k]-[c_k] |^2 を最小にする [A_i] は
  [A_i] = [A_i(0)] + ( t[∂C_k,i] [∂C_k,i] )^(-1) t[∂C_k,i] ( [c_k] - [C_k(0)] )
で求めます。( これを正規方程式というが、次元が大きい場合は正規方程式をそのまま使うと精度が悪くなるので、三角分解した平方根フィルターを使う )
この方が F_i の微分 ( [C_k] の 2 階微分 ) を計算しなくてすむので楽ですが、状況が悪いと解の周りで振動して収束しなくなるので、加速ではなく減速係数を使う必要があったりします。

投稿: hirota | 2007年7月18日 (水) 17時46分

この記事へのコメントは終了しました。

トラックバック


この記事へのトラックバック一覧です: 非線形最小二乗法(JEA式の作成過程):

« 電子を大きさのない点であると考える背景 | トップページ | 安倍政権に思う »