34

私はscipy.optimize.leastsqいくつかのデータを適合させるために使用しています。これらの推定値の信頼区間を取得したいので、cov_x出力を調べますが、これが何であるか、およびこれからパラメーターの共分散行列を取得する方法について、ドキュメントは非常に不明確です。

まず、それはジャコビアンであると書かれていますが、メモには「ヘッセ行列のジャコビアン近似である」と書かれcov_xているため、実際にはジャコビアンではなく、ジャコビアンからの近似を使用したヘッセ行列です。これらのステートメントのどれが正しいですか？

第二に、私にとってこの文は紛らわしいです：

パラメータ推定値の共分散を取得するには、この行列に残差分散を掛ける必要があります。を参照してくださいcurve_fit。

私は確かにcurve_fit彼らがどこで行っているのかソースコードを見に行きます：

s_sq = (func(popt, *args)**2).sum()/(len(ydata)-len(p0))
pcov = pcov * s_sq

これは乗算に対応しcov_xますs_sqが、この方程式はどのリファレンスにも見つかりません。誰かがこの方程式が正しい理由を説明できますか？私の直感ではcov_x、派生物（ジャコビアンまたはヘシアン）であると想定されているので、逆にする必要があると教えてくれたので、私は考えていました： 私が欲しいものはcov_x * covariance(parameters) = sum of errors(residuals)どこにsigma(parameters)ありますか。

どうすればcurve_fitが行っていることを、たとえば、私が見ているものと結び付けることができますか。ウィキペディア： http： //en.wikipedia.org/wiki/Propagation_of_uncertainty#Non-linear_combinations

4

3 に答える 3

29

OK、私は答えを見つけたと思います。最初の解決策：cov_x * s_sqは、必要なパラメーターの共分散です。対角要素の平方根を取ると、標準偏差が得られます（ただし、共分散に注意してください！）。

残差分散=カイ二乗の減少=s_sq= sum [（f（x）-y）^ 2] /（Nn）、ここでNはデータポイントの数、nはフィッティングパラメーターの数です。カイ二乗を減らしました。

私の混乱の理由は、leastsqによって与えられるcov_xは、実際には他の場所でcov（x）と呼ばれるものではなく、縮小されたcov（x）または分数のcov（x）であるためです。他の参考文献のいずれにも表示されない理由は、数値計算には役立つが教科書には関係のない単純な再スケーリングであるためです。

ヘッセ行列とジャコビアンについては、ドキュメントの文言が不十分です。ジャコビアンは最小でゼロであるため、明らかなように、両方の場合で計算されるのはヘッセ行列です。彼らが意味するのは、ヘッセ行列を見つけるためにジャコビアンの近似を使用しているということです。

さらなるメモ。Curve_fitの結果は、実際にはエラーの絶対サイズを考慮しておらず、提供されたシグマの相対サイズのみを考慮しているようです。これは、エラーバーが100万倍変化しても、返されるpcovは変化しないことを意味します。もちろんこれは正しくありませんが、標準的な方法のようです。Matlabは、カーブフィッティングツールボックスを使用するときに同じことを行います。正しい手順は次のとおりです：https ：//en.wikipedia.org/wiki/Linear_least_squares_（mathematics）#Parameter_errors_and_correlation

少なくとも線形最小二乗法では、最適なものが見つかったら、これを行うのはかなり簡単なようです。

于 2013-02-13T15:50:40.477 に答える
7

同様の質問を検索しているときにこの解決策を見つけましたが、HansHarhoffの回答にはわずかな改善しかありません。leastsqからの完全な出力は、infodict ['fvec'] = f（x）-yを含む戻り値infodictを提供します。したがって、縮小されたカイ2乗を計算するには=（上記の表記法で）

s_sq = (infodict['fvec']**2).sum()/ (N-n)

ところで。これを解決するために手間のかかる作業のほとんどを行ってくれたHansHarhoffに感謝します。

于 2013-07-07T03:12:11.127 に答える
3

算数

まず、線形回帰から始めます。多くの統計問題では、変数にいくつかの未知のパラメーターを持ついくつかの基礎となる分布があると想定し、これらのパラメーターを推定します。線形回帰では、従属変数yiが独立変数xijと線形関係にあると仮定します。

y i = xi1β1 + ... + xipβp +σεi 、i = 1、...、n 。

ここで、εiは独立した標準正規分布を持ち、βj 'sはp個の未知のパラメーターであり、σも未知です。これは行列形式で書くことができます：

Y =Xβ+σε、

ここで、Y、β、およびεは列ベクトルです。最良のβを見つけるために、二乗和を最小化します

S =（Y-Xβ）T（Y-Xβ）。

解決策を書き出すだけです。

β^=（X T X）-1XTY。

Yを特定の観測データと見なす場合、β^はその観測の下でのβの推定値です。一方、Yを確率変数と見なすと、推定量β^も確率変数になります。このようにして、β^の共分散が何であるかを見ることができます。

Yには多変量正規分布があり、β^はYの線形変換であるため、β^にも多変量正規分布があります。β^の共分散行列は次のとおりです。

Cov（β^）=（X T X）-1 X T Cov （Y）（（X T X）-1 X T）T =（X T X）-1σ2。

ただし、ここではσが不明であるため、推定する必要があります。させたら

Q =（Y-Xβ^）T（Y-Xβ^）、

Q /σ2はn - p自由度のカイ二乗分布を持っていることが証明できます（さらに、Qはβ^に依存しません）。これにより

σ^ 2 = Q /（n --p）

σ2の偏りのない推定量。したがって、Cov（β^）の最終的な推定量は次のようになります。

（X T X）-1 Q /（n-p）。

SciPy API

curve_fitが最も便利です。2番目の戻り値pcovは、β^の共分散の推定値です。これは、上記の最終結果（X T X）-1 Q /（n --p）です。

ではleastsq、2番目の戻り値cov_xは（X T X）-1です。Sの表現から、X T XはSのヘッセ行列（正確にはヘッセ行列の半分）であることがわかります。そのため、ドキュメントcov_xにはヘッセ行列の逆行列が記載されています。cov_x共分散を取得するには、Q /（n --p）を掛ける必要があります。

非線形回帰

非線形回帰では、yiはパラメーターに非線形に依存します。

y i = f（x i、β1 、 ...、 βp ） + σεi。

βjに関するfの偏導関数を計算できるので、ほぼ線形になります。その場合、計算は基本的に線形回帰と同じですが、最小値を繰り返し近似する必要がある点が異なります。実際には、アルゴリズムは、のデフォルトであるLevenberg–Marquardtアルゴリズムなどのより洗練されたものにすることができますcurve_fit。

シグマの提供に関する詳細

このセクションでは、のとパラメータについて説明sigmaしabsolute_sigmaますcurve_fit。curve_fitYの共分散について事前の知識がない場合の基本的な使用法については、このセクションを無視してかまいません。

アブソリュートシグマ

上記の線形回帰では、y iの分散はσであり、不明です。分散がわかっている場合。パラメータとセットをcurve_fit介してに提供できます。sigmaabsolute_sigma=True

sigma提供された行列がΣであると仮定します。すなわち

Y〜N（Xβ、Σ）。

Yは、平均Xβと共分散Σの多変量正規分布を持ちます。Yの尤度を最大化する必要があります。Yの確率密度関数から、最小化と同等です。

S =（Y-Xβ）TΣ - 1（Y-Xβ）。

解決策は

β ^= （ XTΣ - 1X ）-1XTΣ - 1Y。

と

Cov （β^）=（XTΣ - 1X）-1。

上記のβ^およびCov（β^）は、curve_fitwithの戻り値ですabsolute_sigma=True。

相対シグマ

y iの正確な分散がわからない場合もありますが、異なるy i間の相対的な関係はわかっています。たとえば、y2の分散はy1の分散の4倍です。次に、パスして設定できます。sigmaabsolute_sigma=False

今回

Y〜N（Xβ、Σσ）

既知の行列Σが提供され、未知の数σが与えられます。最小化する目的関数は、σが定数であるため絶対シグマと同じであり、したがって推定量β^は同じです。しかし、共分散

Cov（β^）= （XTΣ - 1X）-1σ2、

未知のσが含まれています。σを推定するには、

Q =（Y-Xβ^）TΣ - 1（Y-Xβ^）。

ここでも、Q /σ2は、n-p自由度のカイ2乗分布を持っています。

Cov（β^）の推定は次のとおりです。

（XTΣ - 1X ）-1 Q /（n-p）。

そして、これはcurve_fitwithの2番目の戻り値ですabsolute_sigma=False。

于 2020-11-14T10:40:26.227 に答える