今日は窒素酸化物濃度(concentrasion of NOX)などを計算する解析的な低煙源拡散式(analytic diffusion formula from low-height source)のパラメータの推定(prediction of parameters)のために用いたことがある非線形最小二乗法(non-linear least square method)について解説してみたいと思います。(直角風のJEA式)
地上の道路上の自動車(automobiles on the road)から発生する窒素酸化物などを,単位長さ当たり(per unit length)の発生源(source)がＱ(Ｎｍ3/(ｍｓ)のｙ軸上で有効煙源高さ(effective height):ｈe＝0(ｍ)の無限線煙源(infinite line-souce)としてモデル化(modelize)し,道路に直角のｘ方向に風が吹いている設定で(wind direction is x normal to road=y-axis),風速(wind speed)も拡散係数(diffusion coefficient)も高さ(height from ground)ｚ(ｍ)のベキ乗則に従って鉛直上方に向かって増加する(increase by power of z)というモデルでの拡散方程式(diffusion equation)の解(solution)を,垂直距離(distance)ｘ(ｍ)と高さｚ(ｍ)の関数と考えて濃度Ｃ(ｘ,ｚ)で表すと,これはＡをＱに比例するあるパラメータ(A is a parameter proportional to Q)として,Ｃ(ｘ,ｚ)＝Ａｘ-sexp(－Ｂｚp/ｘ)という式(Robertsの式)で表わされる,ことがわかります。
いろいろな環境(circumstance)の下で数回にわたって行なわれた実際の実験(experiments)時の煙源長さ(source-length)はもちろん有限(finite)なのですが,とりあえずこれを無限大長さ(infinite length)で近似します。実は誤差関数(Gaussian error function)でもって有限長さの効果(finite-legth effect)を取り入れることもできるのですが,今日の話題では割愛します。実験はトレーサーガス(tracer gas)を地上にある有限長さの直線状のパイプ(linear pipe)に均等間隔(equal distance)で開けた穴(holes)から拡散させるというものです。
そのパイプ状の発生源をｙ軸としたとき,ｎ個の観測点座標(coordinate of n-obserbation points):(ｘi,ｚi)(ｉ＝1,2,．．．ｎ)において,風向が直角に近い環境のときの実測濃度(measured concentration):ｃiと先に設定した低煙源拡散式による計算値(calculated value):Ｃi＝Ｃ(ｘi,ｚi)とを比較することによって,逆にその式のパラメータＡ,ｓ,Ｂ,ｐを推定します。
具体的には,Ｃ(ｘ,ｚ)＝Ａｘ-sexp(－Ｂｚp/ｘ)の右辺をｆ(Ａ,ｓ,Ｂ,ｐ；ｘ,ｚ)と書いて,Ｃi＝ｆ(Ａ,ｓ,Ｂ,ｐ；ｘi,ｚi)とし,計算値と実測値の誤差の二乗和(sum of square of error):Ｓ(Ａ,ｓ,Ｂ,ｐ)＝∑(Ｃi－ｃi)2を最小(minimum)にするパラメータＡ,ｓ,Ｂ,ｐを求めるわけです。これは次に述べる非線形最小二乗法のパラメータ数が４個のときの例(example)になっています。
ここではより一般的にパラメータがｒ個あるとし,それらを順にＡ1, Ａ2,．．Ａr とします。上のケース(case)ではｒ＝4で,Ａ1＝Ａ,Ａ2＝ｓ,Ａ3＝Ｂ,Ａ4＝ｐです。そして,Ｃi＝Ｃ(ｘi,ｚi)＝ｆ(Ａ1,Ａ2,．．Ａr；ｘi,ｚi)とし,Ｓ(Ａ1,Ａ2,．．Ａr)≡∑(Ｃi－ｃi)2と定義するわけです。
誤差の二乗和Ｓが最小となる条件は通常の線形回帰(linear regression)の計算の場合と同様,Ｆi(Ａ1,Ａ2,．．Ａr)≡∂Ｓ/∂Ａi＝2∑(Ｃk－ｃk)∂Ｃk/∂Ａi＝0 (ｉ＝1,2,．．ｒ)で与えられますから,この r 個の連立方程式(simultaneous equations)を解くことにより,未知数(unknown numbers)Ａ1,Ａ2,．．Ａrを求めることが主目的となります。これらの方程式は一般に非線形ですから,こうした非線形回帰によるパラメータ推定の方法を非線形最小二乗法と呼ぶことにします。
具体的には,Ｆi(Ａ1,Ａ2,．．Ａr)＝0 (ｉ＝1,2,．．ｒ)を線形近似(linear approximation)することにより多変数のニュートン法(Newtonian method of many variables)を実行します。
そこで,まず連立方程式を線形近似します。すなわち,初期値(initial values)としてＡj＝Ａj0 を適当に与えた後に, 0 ＝ Ｆi(Ａ1,Ａ2,．．Ａr)＝Ｆi(Ａ10,Ａ20,．．Ａr0)＋∑(∂Ｆi/∂Ａj)|Ａ＝Ａ0(Ａj－Ａj0)と近似します。
これを行列形式(matrix form)で書くため,Ｄという行列(matrix)をＤ＝(ｄij)≡(∂Ｆi/∂Ａj)で定義し,特にＤ0＝(ｄij0)≡(∂Ｆi/∂Ａj)|Ａ＝Ａ0とします。そして,列ベクトルＡをＡ≡t(Ａ1,Ａ2,．．Ａr)で定義し,特にＡ0≡t(Ａ10,Ａ20,．．Ａr0)として,さらにＦi(Ａ10,Ａ20,．．Ａr0)を成分とする同様な列ベクトルをＦ0と定義すれば,先の線形近似は 0＝Ｆ0＋Ｄ0(Ａ－Ａ0)と表わされます。
これを単純に解くと,Ａ＝Ａ0－Ｄ0-1Ｆ0と近似されることになります。ここで,Ｄ0-1はもちろんＤ0の逆行列(inverse matrix)です。
得られたＡ＝Ａ0－Ｄ0-1Ｆ0を,改めてＡ0として代入(substitute)してＤ0-1Ｆ0を計算し,Ａ＝Ａ0－Ｄ0-1Ｆ0を収束(converge)するまで繰り返し(iterate)ます。すなわち,Ａm+1＝Ａm－Ｄm-1Ｆmの漸化式(recurrence formula)において誤差(error)|Ａm+1－Ａm|の相対値(rate)が十分小さく(sufficiently small)なるまで繰り返し計算すればＡm＝t(Ａ1m,Ａ2m,．．Ａrm)のｍ→∞ での極限値(limit)としてＦi(Ａ1,Ａ2,．．Ａr)＝0 (ｉ＝1,2,．．．ｒ)を満足するＡ＝t(Ａ1,Ａ2,．．Ａr)が得られるというわけです。
さらに実際には収束を速くするために加速係数(acceleration coefficient):ｗを与えてＡm+1＝Ａm－ｗＤm-1Ｆmとする方法(method)を用いたほうがいいと思います。
実際のｒ＝4 の計算では,かつてFortranを使ってプログラミング(programming)しましたが,初期値を適切に取ると各環境のケースについて大体十数回の繰り返し計算でパラメータの推定値が得られました。そして濃度の実測値とその推定パラメータによる計算値との相関係数(correlation coefficient)としては,0.9前後の値が得られ,回帰係数(regression coefficient)の値も0.8 から 1.2程度となって,仮定した計算式が良い近似になることがわかりました。
http://fphys.nifty.com/ (ニフティ「物理フォーラム」サブマネージャー）
http://blog.with2.net/link.php?269343 (ブログ・ランキングの投票）↑ここをクリックすると投票したことになります。