【はじめに】

　

　いうまでもなく，数値解析とはノイズによって歪められた情報の中からシグナルを検出し，見えなかったものが見えてくるようにする方法であり，統計解析とは基準値からのずれを本物とみなすべきか統計的誤差範囲内とみなすべきかを客観的な評価基準によって判定する方法である．

　

　現在，主流となっている解析理論は線形・正規現象を取り扱うものであるが，線形・正規的な手法で現実のデータを解析しようとすると，理論と実際のデータが食い違いをみせることがしばしば経験される．これは旧態依然の嫌いがある旧来の解析法の欠陥といえるのだが，線形関数や正規分布が適合しないにもかかわらず，これらの単純な関数を無理やり用いて何らかの結論をデータから導きだそうとすると，データのもつ情報を十分には活用できなくなる危険性をはらんでいるのである．

　

　データを解析して科学的な結論を引き出そうとする場合，グラフ解析や統計解析に頼らざるを得ないため，研究の効率と質はそれらを適切に利用する能力に大きく左右される．情報化時代を迎え，より良いデータ解析の必要性がますます高まってきたが，データ解析の最新技術がもつ有用性は多くの研究者に認識されているとはいいがたい．たとえば，「非線形解析」・「非正規解析」は旧来の解析法の穴を埋めてくれる新しい解析法であると期待されているにもかかわらず，理解されはじめたのはごく最近のことである．

　

　また，測定データには必ず何らかの誤差がつきまとうものだが，データ解析を適切に行えば最小の実験費用（実験にかかる時間とコスト）の下で得られるべき情報量を最大にすることができるようになり，作業能率は向上する．そのためには誤差の扱い方，データ解析の信頼性について詳しく知っておくことが必要になる．

　

　これまでの「非線形解析」・「非正規解析」は，情報のもつ平均的な指標だけを抽出するものであったが，今後は，誤差解析によってその精度を保証すること，すなわち「精度保証つき計算」が要求されているのである．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　コラム「中心極限定理不要論」では，「非正規解析」における「精度保証つき計算」について述べた（つもりである）．測定値は近似的に正規分布にしたがうと仮定されているが，実際の測定結果は必ずしも正規分布にしたがうものではない．しかしながら，母集団が正規分布でないときであっても，中心極限定理により，標本平均値の分布は測定回数が増えるにつれて正規分布に近づく．−−−中心極限定理は正規分布のもつ重要性を物語っていて，このことを突き詰めていったのが「非正規解析」である．

　

　「非正規解析」における精度保証は，母数の分布を正規分布に還元することによって得られたのであるが，「非線形解析」にとっても同じことがいえる．すなわち，非線形式を線形化することによって精度を与えるのである．今回のコラムでは，「非線形解析」における「精度保証つき計算」について解説することにするが，具体的には，非線形モデルを多変量テイラー展開して，１次の項を除く高次項を無視して線形モデル近似を行うことになる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【y=f(x,θ)型回帰】

　

　ある現象の測定データがｎ個のデータ組（ｘ1，ｙ1），（ｘ2，ｙ2），・・・，（ｘn，ｙn）で与えられ，ｘiとｙi（ｉ＝１〜ｎ）の間には何らかの関数関係があるとしよう．２次元データ（ｘi，ｙi）に何らかのモデル式をあてはめるというデータ処理において，通常用いられるモデル式はｍ個の未知のパラメータθ＝（θ1，θ2，θ3，・・・，θm）を含む陽関数ｙ＝ｆ（ｘ，θ）である．以下の議論では，y=f(x,θ)型回帰は縦軸のデータにだけ誤差がある場合で，横軸のほうの誤差はまったくないかあるいは無視できるという前提をおくことにする．

　

　この関数を真の値の近似値θ0のまわりでテイラー展開し，パラメータの推定値がよければ２次以上の導関数を無視できるため，１次までの項をとると，

　

　　ｆ（ｘ，θ）≒ｆ（ｘ，θ0）＋Σ∂ｆ（ｘ，θ0）／∂θk（θk−θk0）

　

　このように線形化さえできればあとはしめたものであり，その際，精度を保証してくれるのが，誤差伝播の法則（propagation of error）である．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【誤差伝播の法則】

　

　簡単のため，ｙが２つの母数θ1，θ2の関数：ｙ＝ｆ（θ1，θ2）のような関数関係があるとする．ｆ（θ1，θ2）を点（θ10，θ20）のまわりでテイラー展開すると，

　　ｙ−ｙ0＝∂ｆ／∂θ1（θ1−θ10）＋∂ｆ／∂θ2（θ2−θ20）

　

　ここで，（ｙ−ｙ0）^2の期待値を計算すると

　　Ｅ（ｙ−ｙ0）^2

　＝Ｅ（∂ｆ／∂θ1（θ1−θ10）＋∂ｆ／∂θ2（θ2−θ20））^2

　＝（∂ｆ／∂θ1）^2Ｅ（θ1−θ10）^2＋（∂ｆ／∂θ2）^2Ｅ（θ2−θ20）^2＋２（∂ｆ／∂θ1）（∂ｆ／∂θ2）Ｅ（θ1−θ10）（θ2−θ20）

となる．

　

　　Ｅ（θ1−θ10）^2＝varθ1＝cov（θ1，θ1）

　　Ｅ（θ2−θ20）^2＝varθ2＝cov（θ2，θ2）

　　Ｅ（θ1−θ10）（θ2−θ20）＝cov（θ1，θ2）

であり，この式は

　

　varｙ＝（∂ｆ／∂θ1）^2varθ1＋（∂ｆ／∂θ2）^2varθ2

　　　　　＋２（∂ｆ／∂θ1）（∂ｆ／∂θ2）cov（θ1，θ2）

　

で与えられるというのが誤差伝播の法則である．このとき，varθ1,varθ2,cov（θ1，θ2）は尤度方程式の分散共分散行列より求めることができる．

　

　ここで注意したいことは，教科書にある誤差伝播の法則では，右辺第３項

　　２（∂ｆ／∂θ1）（∂ｆ／∂θ2）cov（θ1，θ2）

が省かれて記載されていることが多いのだが，母数θ1，θ2間には通常強い相関があり，右辺第３項を省略すると誤差が過大に評価されてしまう．したがって，母数θの誤差Δθを求めるには，右辺第３項を省略してはならない．

　

　また，θ1とθ2の分散をそれぞれ（Δθ1）^2，（Δθ2）^2，共分散をΔθ1Δθ2と略記することにすると，

　　（Δｙ）^2

　＝（∂ｆ／∂θ1）^2（Δθ1）^2＋（∂ｆ／∂θ2）^2（Δθ2）^2

　　　　　＋２（∂ｆ／∂θ1）（∂ｆ／∂θ2）Δθ1Δθ2

　これは，形式的に

　　（Δｙ）^2＝（∂ｆ／∂θ1Δθ1＋∂ｆ／∂θ2Δθ2）^2

と書くことができる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

　一般に，ｙが母数θ1のみの関数，あるいは，三つ以上の母数の関数のときも同様に扱うことができて，誤差伝播の法則は，以下のように書くことができる．

　

　　（Δｙ）^2＝ΣiΣj（∂ｆ／∂θi）（∂ｆ／∂θj）cov（θi，θj）

　

　誤差の伝播法則は近似式であって，正確な不偏推定値とはいえないが，その証明にはテイラー展開しか用いていないことからもわかるように，θ1，θ2が正規母集団からの標本でない場合でも使えて，いろいろな場面に応用できる．もしも，θ1とθ2が厳密に正規分布に従い，その関数ｙ＝ｆ（θ1，θ2）が簡単な関係式で表されるときは誤差の伝播公式よりも正確な誤差の伝播式を求めることは可能である→【参】．しかし，この公式のほうが分布の形にかかわらず成立するので利用価値は高いと考えられよう．

　

【参考文献】吉沢康和「新しい誤差論」共立出版，第７章

　　　　　　粟屋　隆「データ解析」学会出版センター，第７章

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【y=f(x,θ)型回帰の信頼区間】

　

　次に，y=f(x,θ)型回帰曲線の信頼区間であるが，その信頼区間表示については，統計パッケージ「耕太郎」にすでにサポートされていて，食中毒原因菌の増殖曲線から食中毒の発生を予測する場合など，医学分野においてもしばしば利用されている．

　

　話を最も単純な場合に限定するため，１次回帰ｙ＝ａ＋ｂｘにおいてｘがある値をとるとき，回帰直線の両側にｙの推定値についての信頼区間を求めてみることにするが，教科書にある記載では，詳しく書かれている場合であっても，

　

　（１）回帰直線の信頼区間

は回帰直線のまわりのｙの変動を表わす標準誤差σ／√ｎに，回帰係数の推定誤差分σ√（ｘ−ｘm）^2 ／Σ（ｘi −ｘm）^2が加わって

　　σ√（１／ｎ＋（ｘ−ｘm）^2／Σ（ｘi −ｘm）^2）

　（２）観測値の信頼区間

には更に標準偏差σのずれが加わって

　　σ√（１／ｎ＋（ｘ−ｘm）^2／Σ（ｘi −ｘm）^2 ＋１）

に対応した値をとる．

　（３）明らかに，前者に比べ後者の信頼区間の幅は広くなり，どちらもｘがｘに等しいとき区間の幅がもっとも狭く，ｘから遠ざかるにつれて広くなる双曲線になる．

　（４）直線全体に対する信頼区間はこの検定量が自由度ｎ−２のｔ分布に従うことを利用して求めることができる．

　

と，いかにもそっけない．これでは何のことだかわからないし，理解できないのも当然であろう．そこで，係数ａ，ｂの分散・共分散を求めてみることにしよう．

　

　計算は省略するが，

　　varａ＝σ^2｛１／ｎ＋ｘm^2／Σ（ｘi −ｘm）^2｝

　　varｂ＝σ^2／Σ（ｘi −ｘm）^2

　　cov（ａ，ｂ）＝−σ^2／Σ（ｘi −ｘm）^2

となる．ここで，誤差伝播の法則を用いると，

　　（Δｙ）^2＝（Δａ＋ｂｘ）^2

　　　　　　　＝（Δａ）^2＋２ΔａΔｂｘ＋（Δｂ）^2ｘ^2

　　　　　　　＝σ^2｛１／ｎ＋（ｘ−ｘm）^2／Σ（ｘi −ｘm）^2｝

となって，（１）と同じ結果が得られる．

　

　非線形回帰式y=f(x,θ)の場合も，このことを敷衍していけばよいのであるが，　　E[x]=μ

　　V[x]=E[x^2]-E[x]^2=σ^2

であるとき，確率変数xの分布に関わらず，任意の定数a,bに対して

　　E[ax+b]=aE[x]+b

　　V[ax+b]=a^2V[x]

が成り立つこと，また，回帰では変数と定数の関係が逆転していることを踏まえて，以下のようにステップアップ式に考えてみよう．

　

【問】線形式：ｙ＝ａ1ｆ1（ｘ）＋ａ2ｆ2（ｘ）＋・・・＝Σａkｆk（ｘ）

ただし，ｆk（ｘ）はｘ^2とかｓｉｎｘといったｘの関数であればなんでもよいのであるが，未知パラメータが含まれていてはいけない．この場合，ｙの誤差Δｙは？

（答）（Δｙ）^2＝ΣΣｆi（ｘ）ｆj（ｘ）cov（ａi，ａj）

　　　　　　　　＝ΣΣ（∂ｙ／∂θi）（∂ｙ／∂θj）cov（θi，θj）

　

【問】擬非線形式：ｇ（ｙ）＝Σａkｆk（ｘ）．この場合のｙの誤差Δｙは？

（答）Ｙ＝ｇ（ｙ）＝Σａkｆk（ｘ）とおくと，

　　（ΔＹ）^2＝ΣΣｆi（ｘ）ｆj（ｘ）cov（ａi，ａj）

また，誤差伝播の法則より

　　（ΔＹ）^2＝（∂ｇ／∂ｙ）^2（Δｙ）^2

これより，Δｙが求まる．

　

【問】非線形式：ｙ＝ｆ（ｘ，θ1，θ2，θ3，・・・，θm）において，ｙの誤差Δｙは？

（答）非線形，すなわち，未知の係数がべき乗の形になってたり，指数関数や三角関数の中に含まれていたりする関係があるときでも，非線形式は

　　ｙ−ｙ0≒Σ∂ｆ／∂θkΔθk

と線形化できるから，

　　（Δｙ）^2＝ΣiΣj（∂ｆ／∂θi）（∂ｆ／∂θj）cov（θi，θj）

となることがわかるのである．

　

　この式は線形の場合に一致することに注意されたい．もちろん，線形ではないので厳密には正しくはないが，違いは大きいものではないから，それほどこだわる必要はないだろう．なお，∂ｆ／∂θi，∂ｆ／∂θjは，ともにｘの関数である．また，関数y=f(x,θ)の100(1-α)%の信頼区間は±t(α/2,n-m)Δｙで表されることになる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【最小２乗法の拡張】

　

　２次元データ（ｘi，ｙi）に陽関数ｙ＝ｆ（ｘ）をあてはめるというデータ処理において，ｘの誤差が無視できる場合，残差２乗和ｓ＝Σ（ｙi−ｆ（ｘi））^2，あるいは，これに重みをつけたｓ＝Σｗi（ｙi−ｆ（ｘi））^2を最小にする関数近似が最小２乗法である．

　

　非線形最小２乗法の教科書はこれまで多数出版されているが，大抵はモデル式が非線形陽関数の場合のみを考えて解説している．最小２乗法はさまざまな方法で一般化できるのだが，この目的関数でよしとして済ませてしまうことができるのはもっとも簡単な場合に限られている．つまり，一口に最小２乗法といっても，データが得られた状況によってそれぞれに最適な手法が用意されていて，それらを適宜選択し使い分ける必要が出てくるのである．

　

　最小２乗法の解析条件をおおよそ列挙してみると，

ａ）モデル式は線形か非線形か

　　モデル式は陽関数，陰関数，微分方程式のいずれか

ｂ）横軸側のデータに誤差はあるのかないのか

ｃ）重みつき最小２乗法かどうか

　　重みの中に未知パラメータが含まれるか

ｄ）パラメータの制約条件の有無はどうか

などがあげられる．

　

　本稿ではそのうちの重要なものを取り上げるが，以下，f(x,y)=0型回帰曲線やdy/dx=f(x,y)型回帰曲線のあてはめとその信頼区間を中心として解説を加えることにする．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【f(x,y,θ)=0型回帰】

　

　y=f(x,θ)型回帰は「縦軸のデータにだけ誤差がある場合で，横軸のほうの誤差はまったくないかあるいは無視できる」という前提をおいたものである．その幾何学的な意味は，「実測値と曲線の鉛直方向のユークリッド距離の２乗和が最小になるように未知のパラメータの値を定め，あてはめの誤差を少なくする．」と解釈される．

　

　実用上はこれで十分な場合も多いが，横軸側のデータといえども必ずしも理想的には管理されず，誤差の混入が常であるから，どうしてもｘ，ｙ双方に誤差を考えないといけない場合がある．

　

　たとえば，ある物質と受容体の結合親和性や受容体量を求めるためのラングミュア・プロット（Langmuir plot）では，縦軸・横軸とも測定濃度をプロットするので，ｘ軸側，ｙ軸側ともに誤差があり，その値が大きいほどその誤差の絶対値は大きくなる．また，最高血圧（東京の最高気温）をｘ，最小血圧（仙台の最高気温）をｙとした散布図ではｘがｙの原因である（ｙがｘの原因である）という関係が想定しにくく，通常の回帰の観念を用いるのは適当ではない．

　

　ｘ，ｙ双方に同じ程度の誤差が存在する場合は，測定点からあてはめるべき曲線に下ろした垂線の長さが誤差となり，垂線の長さの平方和が最小になるように定めた曲線がｘとｙの間に成立する量的関係のもっとも確からしい表現と見なすことができる．このような場合，ｘ側もｙ側も対等だから，あてはめる関数を陰関数：f(x,y,θ)=0としておいたほうが都合がよいばかりでなく，円とか楕円とか，多価関数のあてはめも可能になるという利点もある．

　

　横軸側にも誤差がある場合の関数フィッティング法に，粟屋（青山学院大学・理工学部・物理）の２次元フィッティング法がある→【参】．粟屋の関数フィッティングの方法の詳細は記さないが，最尤法に基づいて，垂線の長さ（正確にはマハラノビス距離）の平方和の最小化を一般化したものである．陽関数におけるガウス・ニュートン法を陰関数まで取り扱いを拡張したもので，従来の最小２乗法，すなわち，ｙのみに誤差がある，あるいは，ｘのみに誤差があるときの取り扱いも可能になっていて，最小２乗法の決定版と考えることができる．

　

　なお，ｆ（ｘ，ｙ）＝０を一般的に取り扱うためには，陰関数の描画ルーチンなどが必要となる．その点，f(x,y,θ)=0型回帰は原理的に煩わしいといえるだろう．→【補】

　

【参考文献】粟屋隆「データ解析」学会出版センター，第１２章

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【f(x,y,θ)=0型回帰の信頼区間】

　

　以下，∂f/∂θk，∂f/∂x，∂f/∂yをそれぞれfk,fx,fyと略記するが，この場合も，未知母数の信頼区間を与えてくれる分散共分散行列{cov(θi,θj)}は正規方程式の係数行列の逆行列として求めることができるし，関数ｆの誤差は，パラメータの回りでテイラー展開して，

　　(Δf)2=ΣΣfifjcov(θi,θj)

で与えられることも，y=f(x,θ)=0型回帰の場合と同様である．

　

　横軸側にも誤差がある場合の信頼区間表示で異なっている点は，この誤差を縦軸方向と横軸方向に分解する必要があることである．再び，誤差伝播の法則

　　(Δf)^2=(fxΔx+fyΔy)^2

を用いるが，この式は，いわば誤差の分散公式

　　σ2(x+y)=σ2(x)+σ2(y)+2rσ(x)σ(y)

であって，z=x+yとするとｚ軸方向の合成分散をｘ軸，ｙ軸の２つの方向に分配すると考えることができる．

　

　その際，ベクトル（Δｘ，Δｙ）は法線方向を向くが，陰関数ｆ（ｘ，ｙ）＝０上の点（ｘ，ｙ）で接線の方程式を求めるには，２変数関数の微分の知識が必要で，ｆ（ｘ，ｙ）＝０のｙをｘの関数（ｆ（ｘ，ｙ（ｘ））＝０）とみなして，両辺をｘで偏微分すれば２変数関数の合成微分の公式によって

　　∂ｆ／∂ｘｄｘ／ｄｘ＋∂ｆ／∂ｙｄｙ／ｄｘ＝０

すなわち，ｆx＋ｆyｄｙ／ｄｘ＝０より，

　　ｙ’＝ｄｙ／ｄｘ＝−ｆx／ｆy

が得られる．これが陰関数定理であるが，これにより，

　　Δｙ／Δｘ＝ｆy／ｆx

となる．→【補】

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【dy/dx=f(x,y,θ)型回帰】

　

　たとえば，ｘをある生態系におけるウサギの個体数，これをエサとするキツネの個体数をｙとする．それらの変化は被食者と捕食者の生存闘争モデルとして有名なロトカ・ヴォルテラ（Lotka-Volterra）の微分方程式

　　ｄｘ／ｄｔ＝ａｘ−ｃｘｙ

　　ｄｙ／ｄｔ＝−ｂｙ＋ｃｘｙ　　　（ａ，ｂ，ｃ＞０）

で表される．

　

　この微分方程式中には非線形項ｘｙが含まれ，もはや解析的に解くことはできない．そこで，微分方程式を数値的に解く，すなわち，常微分方程式の初期値問題のパラメータを計算する部分に，数値積分法であるルンゲ・クッタ法などを組み込み，数値的に解きつつ非線形最小２乗法にかけるというやり方で，パラメータの計算を行なうことになる．→【補】

　

　この例はモデル式が１階常微分方程式で与えられる場合のあてはめ例であるが，高階常微分方程式の場合も，適当な変数を介して，多元連立常微分方程式に変形できるため，高階になっても連立１階常微分方程式と同様に解くことができる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【dy/dx=f(x,y,θ)型回帰の信頼区間】

　

　一般に，自然現象や社会現象の数理モデルは１階・高階微分方程式で表わされることが多い．しかし，微分方程式を解析的に解くことができない場合は，微分方程式を数値的に解きつつ関数フィッティングを行うので，その計算量は極めて膨大になることを覚悟せねばならない．

　

　信頼区間に関しては，y=f(x,θ)型やf(x,y,θ)=0型と同様に，微分係数ｄｙ／ｄｘの誤差は

　　(Δdy/dx)^2=ΣΣfifjcov(θi,θj)

で与えられ，この誤差を数値積分法に組み入れることによって信頼区間を計算することができる．その場合，横軸側には誤差はないという前提をおいているので，誤差を縦軸方向にのみ分配すればよい．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【考察】

　

　従来の関数フィッティングでは，横軸のほうの誤差は無視できるとか，モデル関数式が陽に表現されていなければならないという制約があるが，ここに掲げたような方法でこの制約はずいぶんと解消され，モデル関数式をかなり拡張させることができた．

　

　また，関数フィッティングは

　　y=f(x)→f(x,y)=0→dy/dx=f(x,y)

の順にむずかしくなるが，いずれの場合にも，信頼区間を求めるのに，誤差伝播の法則が共通して用いられていることもおわかり頂けたかと思う．

　

　特に，関数フィッティングの応用を考えると，あてはまりの良さを追求するよりも，信頼区間を推定できることのほうが重要なので，ここで示した信頼区間表示法の意義は大きいと考えられる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【間違いだらけのソフト選び】

　

　非線形最小２乗法に関しては，決して多いとはいえないながらも，これまでにいくつかのソフトウェアが開発されてきた．それらはサブルーチンであったり，商用のプログラムパッケージであったりするわけであるが，非線形最小２乗法ソフトの優劣の見分け方について解説しておきたい．

　

　ａ）演算アルゴリズムが充実していること．

　これは解の収束安定性のためには必須であり，少なくとも，ガウス・ニュートン法，マルカール法，シンプレックス法がサポートされていることが最低条件である．現在，ほとんどのソフトはこの点をクリアしている．→【参】

　

　ｂ）パラメータの誤差評価が精確になされていること．

　誤差を求める場合，分散・共分散行列から求めるのだが，ニュートン法が最もよい方法である．ガウス・ニュートン法で計算しても，ニュートン法とさしたる違いはないので，それほどこだわる必要はないだろう．

　しかし，たとえば，ガウス・ニュートン法で求めた誤差とマルカール法で求めた誤差が食い違ったら，そのソフトの誤差評価は不精確だと思ってよい目安となる．また，シンプレックス法で計算した場合に，誤差が表示されればそのソフトの信頼性は高いと考えてよいだろう．→【参】

　

　ｃ）グラフ機能，とくに回帰曲線の信頼区間表示が可能であること．

　パラメータの誤差は求められていても，信頼区間表示機能が付いていないソフトは精度保証が十分なされているとはいえない．

　非線形最小２乗法によるy=f(x)型回帰とその信頼区間表示については，食中毒原因菌の増殖曲線から食中毒の発生を予測する場合など，医学分野においてもしばしば利用されている．したがって，応用上とくにこの点は重要であるが，残念ながらこれをクリアしているソフトはほとんど見あたらない．

　

　どんなコンピュータ・ソフトでも，使い方の正否で誤りが起こるが，非線形最小２乗法プログラムの誤差評価は精確になされていないことも多いので，とくに注意が必要である．プログラムを作るのも人間である以上，コンピュータは誤りを犯すし，問題なのはコンピュータの出力をいつの場合も鵜呑みにする人間である．コンピュータプログラムには間違ったものがいくらでもあるといったら驚かれるかもしれないが，実際，世界中の多くの非線形最小２乗法プログラムは潜在的なバグをもっていて，いまだ顕在化していないのはプログラマー自身も（無論ユーザーも）それに気づいていないだけの話なのである．プログラムのバグに対しても，柔軟で実際的な姿勢が求められる所以である．

　

　また，少なくとも科学を志す者にとっては，たとえ既成のソフトウェアを利用するにしても，どのような処理がなされているのかその理論的背景を知っておくべきで，The computer told me!と盲信する姿勢は，The tail wags a dog!（主客転倒）であり，厳に慎みたいと思う．

　

【参考文献】佐藤郁郎「最小２乗法その理論と実際」山海堂

　　　　　　佐藤郁郎「最小２乗法ソフト耕太郎のすべて」山海堂

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【おわりに】

　

　数値解析・グラフ解析と統計解析の基礎的な方法を使うだけでもさまざまな現象について深い理解が可能になるが，知的基盤をさらに拡大するにはより巧妙な数学的技法を用いて一連の測定結果から有益な情報を最大限取り出す必要がある．これがdatalogy（データ解析学）の目的である．

　

　私は時代に竿をさした風変わりな人間であるから「非線形解析」やら「非正規解析」やら浮き世離れした研究に携わることになったのだが，このようなdatalogyの方法論の確立に加え，コンピュータを用いる科学計算の研究に携わってきて，ここ１０年の間に開発したプログラムは多数におよぶ．本稿で述べた方法論は非線形解析用ソフト「耕太郎」，非正規解析用ソフト「麦」などのパッケージにすでにまとめられている．世に類似のソフトは少なからず出ているとは思うが，私自身がシニカルでエキセントリックな人間だからこそ「耕太郎」も「麦」も少し趣の異なるソフトに仕上がったと自負している．

　

　精度保証つき非線形解析など，ここで解説した方法が提案されたのはそれほど新しくないが，多くの研究者に理解されはじめたのはごく最近のことである．この機会に是非とも理解を深められたい．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】陰関数定理

　

　　ｙ’＝ｄｙ／ｄｘ＝−ｆx／ｆy

が陰関数定理と呼ばれるものであるが，さらにｆx＋ｆyｙ’＝０の式をｘの関数とみて，この両辺をｘで微分すれば

　　ｆxx＋ｆxyｙ’＋（ｆyx＋ｆyyｙ’）ｙ’＋ｆyｙ”＝０

より

　　ｙ”＝ｄ2ｙ／ｄｘ2

　　　　＝−１／ｆy（ｆxx−２ｆxyｆx／ｆy＋ｆyyｆx^2／ｆy^2）

が得られる．決して，ｙ”＝−ｆxx／ｆyyなどというでたらめを書かないように！

　

【問】ｙ＋ｓｉｎｙ＋ｘ^2−ｅ^x＋１＝０において，ｄｙ／ｄｘを求めよ．

（答）ｆx＝２ｘ−ｅ^x，ｆy＝１＋ｃｏｓｙ．したがって，ｄｙ／ｄｘ＝−（２ｘ−ｅ^x）／（１＋ｃｏｓｙ）．

　

【問】ｙ＝ｘ^xとｙ＝ｘ^1/xを微分せよ．

（答）それぞれｙ’＝（ｌｏｇｘ＋１）ｘ^xとｙ’＝（１−ｌｏｇｘ）ｘ^1/x-2．ｘ^1/xはｘ＝ｅのとき最大値１．４４４６・・・をとる．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】微分方程式の解法

　

　物理や化学に限らず，日常の現象は微分方程式で表される場合が多いのですが，微分方程式を解析的に解くことのできる場合は非常に少なく，数値計算により解く必要がでてきます．

　

（１）オイラー法

　取り扱う微分方程式をｄｙ／ｄｘ＝ｆ（ｘ，ｙ）と表現しましょう．ある点ｘ0におけるｙの値が既知でｙ0のとき，

　　ｄｙ／ｄｘ＝ｆ（ｘ0，ｙ0）

これが，Δｘだけｘを変化させていっても成立するものとするとΔｙ／Δｘ＝ｆ（ｘ0，ｙ0）ですから，

　　Δｙ＝Δｘ・ｆ（ｘ0，ｙ0）

　したがって，ｘ0からΔｘ変化した場合のｙの値は

　　ｙ＝ｙ0＋Δｘ・ｆ（ｘ0，ｙ0）

により求まります．ｘをΔｘ分だけ変化させてこの式を繰り返し計算することにより，任意のｘの範囲まで関数ｙの値を求めることができるようになります．

　

　この特殊解を得る方法をオイラー法といいます．すなわち，オイラー法では初期値（ｘn，ｙn）が与えられたとき（ｘn+1，ｙn+1）は

　　ｘn+1＝ｘn＋ｈ

　　ｙn+1＝ｙn＋ｈ・ｆ（ｘn，ｙn）

で求められます．ただし，Δｘをｈで書き換えています．

　

（２）ルンゲ・クッタ法

　オイラー法ではｈの区間を直線で近似していますから，ｈを余程小さくとらない限り誤差が累積してしまうのが欠点です．この欠点を除いた格段に精度の良い方法がルンゲ・クッタ法（４次法・１／６公式）です．その微分方程式の数値解ｙは次式にて表されます．

　　ｙn+1＝ｙn＋（ｋ0＋２ｋ1＋２ｋ2＋ｋ3）／６

ここに，

　　ｋ0＝ｈ・ｆ（ｘn，ｙn）

　　ｋ1＝ｈ・ｆ（ｘn＋ｈ／２，ｙn＋ｋ0／２）

　　ｋ2＝ｈ・ｆ（ｘn＋ｈ／２，ｙn＋ｋ1／２）

　　ｋ3＝ｈ・ｆ（ｘn＋ｈ，ｙn＋ｋ2）

　

　微分方程式は解析的な形に書き下せないことが多く，実際の場面では数値的に解いたほうが見通しがよい場合が少なくありません．オイラー法，ルンゲ・クッタ法の近似法は連続変数を離散変数で置き換えることからはじめ，テイラー展開を有限で打ち切りそれを数値的に積分して解を求めるものです．オイラー法は原理が簡単なのですが，誤差が大きいので実際の数値計算にはほとんど用いられていません．実際のプログラムではルンゲ・クッタ法がもっとも多く用いられる差分解法になっています．

　

　微分方程式の数値解法の際，ｈの値にはとくに注意する必要があります．これが大きすぎると誤差は大きくなり，小さすぎるとまるめ誤差を生じることになるからです．最適のｈは経験的に決められますが，固定きざみではなく，変化量の大きい部分ではｈを小さくし，変化量の小さい部分では大きくするといったように，その都度の変化率によってｈを自動的に修正して計算させる可変きざみも必要になってきます．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝

　

【補】陰関数を描く

　

　ここでは，微分方程式の解法の応用として，一般的な陰関数ｆ（ｘ，ｙ）＝０の描き方について考えてみます．

　

　陽関数の場合，媒介変数ｔを導入してｘ＝ｔ，ｙ＝ｇ（ｔ）とすれば，これはｙ＝ｇ（ｘ）と同じことになるので，陽関数は必ずパラメータ表示に書き直すことができます．これと同様に陰関数ｆ（ｘ，ｙ）＝０において，一意に定まるパラメータｔを導入してｘ＝ｘ（ｔ），ｙ＝ｙ（ｔ）と表わされるものとすると，座標（ｘ，ｙ）の代わりに（ｔ，ｘ（ｔ）），（ｔ，ｙ（ｔ））を用いることができるため，陽関数と同様の取り扱いが可能になります．

　

　すなわち，陽関数，陰関数いずれの場合であってもｘ−ｙ平面上の点（ｘ，ｙ）を第３の変数ｔをパラメータとする点（ｘ（ｔ），ｙ（ｔ））ととらえると，ｔの変化とともに描かれる点の軌跡がグラフになります．このように陰関数でなしに，陽な形あるいは媒介変数表示の形にして書くことができれば，多価関数が自然に定義できる利点があり申しぶんありません．

　

　したがって，陰関数表示された曲線ｆ（ｘ，ｙ）＝０のグラフを描くことはｆ（ｘ（ｔ），ｙ（ｔ））＝０を満たす２つの関数ｘ（ｔ），ｙ（ｔ）を見つけだすことと同等の問題になります．なんとかして媒介変数表示したいものですが，そのためには微分方程式の差分解法の知識が必須なものになってきます．

　

　そこでまず，ｆ（ｘ0，ｙ0）＝０を満たす２つの数ｘ0，ｙ0が求まっているとします．ｘ（０）＝ｘ0，ｙ（０）＝ｙ0とすると点（ｘ0，ｙ0）の近傍の点（ｘ，ｙ）＝（ｘ（Δｔ），ｙ（Δｔ））で，小さなΔｔに対して

　　ｆ（ｘ，ｙ）＝ｆ（ｘ0，ｙ0）＋Δｘｆx（ｘ0，ｙ0）＋Δｙｆy（ｘ0，ｙ0）＋高次の項

が成立します．

　

　ｘ＝ｘ0＋Δｘ，ｙ＝ｙ0＋Δｙ，ｆ（ｘ0，ｙ0）＝０，ｆ（ｘ，ｙ）＝０より

　　Δｘｆx（ｘ0，ｙ0）＋Δｙｆy（ｘ0，ｙ0）≒０

したがって，一般性を失なうことなく，全微分型方程式

　　ｄｆ／ｄｔ＝ｆx（ｘ，ｙ）ｄｘ／ｄｔ＋ｆy（ｘ，ｙ）ｄｙ／ｄｔ＝０

が得られます．

　

　このとき１階の連立微分方程式があって

　　ｄｘ／ｄｔ＝−ｆy（ｘ，ｙ），

　　ｄｙ／ｄｔ＝　ｆx（ｘ，ｙ）

を満たすものとします（ハミルトン型方程式）．ハミルトン型の方程式では

　　ｄｆ／ｄｔ＝ｆx（ｘ，ｙ）ｄｘ／ｄｔ＋ｆy（ｘ，ｙ）ｄｙ／ｄｔ＝０

が成り立ちますから，ｔについて積分するとｆ（ｘ，ｙ）＝ｃ，初期条件ｆ（ｘ0，ｙ0）＝０よりｃ＝０，すなわちｆ（ｘ，ｙ）＝０となります．この式の場合，ｆはｔに関係なく一定で，点（ｘ，ｙ）は等高線上を進むと考えられますから，軌道に沿ってエネルギーが変化しない系（保存系）を表わすものと解釈されます．このような関数ｆをハミルトン関数といい，物理の世界では運動の全エネルギーを表わすものとして有名です．

　

　以上により，陰関数のグラフを描くことはハミルトン型の２元連立常微分方程式

　　ｄｘ／ｄｔ＝−ｆy（ｘ，ｙ），

　　ｄｙ／ｄｔ＝　ｆx（ｘ，ｙ）

および

　　ｄｘ／ｄｔ＝　ｆy（ｘ，ｙ），

　　ｄｙ／ｄｔ＝−ｆx（ｘ，ｙ）

を初期条件ｘ（０）＝ｘ0，ｙ（０）＝ｙ0のもとで解くことに帰着されますが，前出のルンゲ・クッタ法による差分解法を用いると，この解となる軌道を求めることができます．

　

　最後に，ｆ（ｘ0，ｙ0）＝０を満たす２つの数ｘ0，ｙ0を求めてみます．この解はｘ0を適当にとりｆ（ｘ0，ｙ）をｙの関数と考えて，ニュートン・ラプソン法を用いて解ｙ＝ｙ0を求めることができます．２元連立方程式

　　ｆ（ｘ，ｙ）＝０，

　　ｙ＝ｙ0

においては，更新公式Δｘ＝−ｆ（ｘ，ｙ）／ｆx（ｘ，ｙ），Δｙ＝０が得られます．この式は１変数関数におけるニュートン法の公式

　　Δｘ＝−ｆ（ｘ）／ｆ’（ｘ）

において微分を偏微分に変えただけで全く同じ形式をとっています．

　

＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝＝