マトリックス

数学はNonLinearModelFitでCovarianceMatrixをどのように計算しますか



How Does Mathematica Calculate Covariancematrix Nonlinearmodelfit



解決：

注：私が生まれるずっと前から、線形モデルの計画行列と呼ばれていました。統計学のクラスは、用語に役立ちます。

ただし、非線形設定では、モデルが線形の場合でも通常の計画行列になる、少し複雑な項目であることがわかります。



重み付けの問題は無視します。 NS 非線形モデル用です。以下に、を使用した非線形モデルの例を示します。NonlinearModelFitと、基礎となる計算を表示するブルートフォース最尤法。

（*非線形モデルからデータを生成します*）SeedRandom [12472475]; x = Flatten [Table [{0、1、2、3、4、5、6、7、8}、{i、5}]]; n =長さ[x]; aa = 1; bb = 0.5; σσ= 0.05; y = aa Exp [-bb x] + RandomVariate [NormalDistribution [0、σσ]、n]; （* NonlinearModelFitに適合*）nlm = NonlinearModelFit [Transpose [{x、y}]、Exp [-b z]、{a、b}、z]; （*ブルートフォース法に適合*）（*尤度の対数を取得*）logL = LogLikelihood [NormalDistribution [0、σ]、y --a Exp [-b x]]; （*パラメーターの最尤推定値を検索します*）mle = FindMaximum [{logL、σ> 0}、{{a、1}、{b、0.5}、{σ、0.05}}]; （*ここで、共分散行列を近似してから推定します*）（* Xの定義、計画行列に注意してください*）X = D [a Exp [-b x]、{{a、b}}]; （* n /（n-2）による乗算は、σ^ 2の不偏推定を使用することです*）（* n-2は、実際にはサンプルサイズからパラメーターの数を引いたものです*）cov =（n /（n-2））（ σ^ 2Inverse [Transpose [X]。X]）/。 mle [[2]]; （*パラメーター推定値*）nlm ['BestFitParameters']（* {a-> 0.9842981255137611、b-> 0.5023511174255275} *）mle [[2]] [[1 ;; 2]]（* {a-> 0.9842981255018496、b-> 0.50235111740454} *）（*共分散行列の推定値*）nlm ['CovarianceMatrix']（* {{0.00034959357568038475,0.00016592363407431294}、{0.00016592363407431294,0.0002923502562857745}} *）cov （* {{0.00034959357700581665,0.0001659236346998225}、{0.0001659236346998225,0.00029235025737224096}} *）Rationalize [cov / nlm ['CovarianceMatrix']、0.00001]（*σ^ 2の推定値*）nlm ['EstimatedVariance']（* 0.002014630678417603 *） （（n /（n-2））σ^ 2 /。mle[[2、3]]）（* 0.0020146306860807216 *） 
途中で答えます。だから私は計画行列を誤解しました、それはDesignMatrix OF X. ok、Xと呼ばれるのでばかげています。計画行列の形式は{{1.、Xestimate_1}、{1.、Xestimate_2}、...}です。次のコードは、と同じ行列を生成しますnlm ['CovarianceMatrix']はそうします。ただし、線形モデルの場合のみb + a * x。したがって、まだ未解決の問題があります。 :-)
 
 
 
f =関数[{a、b、x}、b + a * x];データ= {{6.、21。}、{8.、31。}、{9.、39。}};重み= {0.5、4.5、0.1}; nlm = NonlinearModelFit [data、f [a、b、x]、{a、b}、x、Weights-> weights、VarianceEstimatorFunction-> Automatic] yestimates = Table [nlm [data [[i、1]]]、{ i、Length [data]}] xestimates = Flatten [Table [x /。 Solve [nlm [x] == yestimates [[i]]]、{i、Length [data]}]] xdesignmatrix = Table [{1.、xestimates [[i]]}、{i、Length [data]} ]; Varianceestimate = nlm ['EstimatedVariance']; matrixofweights = DiagonalMatrix [weights]; Varianceestimate * Inverse [Transpose [xdesignmatrix] .matrixofweights.xdesignmatrix]％// MatrixForm nlm ['CovarianceMatrix'] // MatrixForm 
 
出力：{{25.9587、-3.29752}、{-3.29752、0.421488}}