scipyでカーブフィッティングの可能性を計算するにはどうすればよいですか?

StackOverflow https://stackoverflow.com//questions/23004374

質問

次のような非線形モデルの適合があります。

Curve fit

濃い実線はモデルの適合、灰色の部分は生データです。

質問の短いバージョン:対数尤度比検定を実行できるように、このモデルが適合する可能性を取得するにはどうすればよいですか?残差は正規分布すると仮定します。

私は統計については比較的初心者ですが、現在の考えは次のとおりです。

  1. カーブフィットから残差を取得し、残差の分散を計算します。

  2. この方程式を使用してくださいLog-likelihood for normal distributionsそして、残差の分散をシグマ二乗に代入し、x_i を実験として、mu をモデルの適合として使用します。

  3. 対数尤度比を計算します。

これら 2 つのフルバージョンの質問について、誰か助けてくれませんか?

  1. 私のやり方は正しいのでしょうか？（そう思いますが、確認していただけると本当に助かります！）

  2. python/scipy/statsmodels にこれを行うための既製の関数はありますか?

役に立ちましたか？

解決

あなたの尤度関数

enter image description here

これは、単にガウス分布の確率密度関数の対数の合計です。

enter image description here

の可能性があります 剰余にミューとシグマを当てはめる, の可能性ではありません データを考慮したモデル. 。あなたのアプローチを一言で言えば、 間違っている.

@usethedeathstarがすでに述べたことに従って、非線形最小二乗を実行しているため、まっすぐに進む必要があります F-test. 。。以下を変更した次の例を考えてみましょう。 http://www.walkingrandomly.com/?p=5254, 、そして私たちは実施します F-test を使用して R. 。そしてそれをどのように翻訳するかについて議論します python 最後に。

# construct the data vectors using c()
> xdata = c(-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9)
> ydata = c(0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001)
# some starting values
> p1 = 1
> p2 = 0.2
> p3 = 0.01

# do the fit
> fit1 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata), start=list(p1=p1,p2=p2))
> fit2 = nls(ydata ~ p1*cos(p2*xdata) + p2*sin(p1*xdata)+p3*xdata, start=list(p1=p1,p2=p2,p3=p3))

# summarise
> summary(fit1)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1 1.881851   0.027430   68.61 2.27e-12 ***
p2 0.700230   0.009153   76.51 9.50e-13 ***
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08202 on 8 degrees of freedom

Number of iterations to convergence: 7 
Achieved convergence tolerance: 2.189e-06

> summary(fit2)

Formula: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata

Parameters:
   Estimate Std. Error t value Pr(>|t|)    
p1  1.90108    0.03520  54.002 1.96e-10 ***
p2  0.70657    0.01167  60.528 8.82e-11 ***
p3  0.02029    0.02166   0.937     0.38    
---
Signif. codes:  0 ?**?0.001 ?*?0.01 ??0.05 ??0.1 ??1

Residual standard error: 0.08243 on 7 degrees of freedom

Number of iterations to convergence: 9 
Achieved convergence tolerance: 2.476e-06

> anova(fit2, fit1)
Analysis of Variance Table

Model 1: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata) + p3 * xdata
Model 2: ydata ~ p1 * cos(p2 * xdata) + p2 * sin(p1 * xdata)
  Res.Df Res.Sum Sq Df     Sum Sq F value Pr(>F)
1      7   0.047565                             
2      8   0.053813 -1 -0.0062473  0.9194 0.3696

ここには 2 つのモデルがあります。 fit1 には 2 つのパラメーターがあるため、剰余には 8 つの自由度があります。 fit2 には追加パラメータが 1 つあり、剰余には 7 自由度があります。モデル 2 は大幅に優れていますか?いいえ、F 値は 0.9194 です。 (1,7) 自由度は大きくありませんが、それは重要ではありません。

ANOVA テーブルを取得するには:残留DFは簡単。剰余平方和: 0.08202*0.08202*8=0.05381 そして 0.08243*0.08243*7=0.04756293 （知らせ： '残存標準誤差:7 自由度では 0.08243', 、など）。で python, 、次の方法で入手できます。 (y_observed-y_fitted)**2, 、 以来 scipy.optimize.curve_fit() 残差は返しません。

の F-ratio は 0.0062473/0.047565*7 そして P 値を取得するには: 1-scipy.stats.f.cdf(0.9194, 1, 7).

それらを組み合わせてください python 同等：

In [1]:

import scipy.optimize as so
import scipy.stats as ss
xdata = np.array([-2,-1.64,-1.33,-0.7,0,0.45,1.2,1.64,2.32,2.9])
ydata = np.array([0.699369,0.700462,0.695354,1.03905,1.97389,2.41143,1.91091,0.919576,-0.730975,-1.42001])
def model0(x,p1,p2):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)
def model1(x,p1,p2,p3):
    return p1*np.cos(p2*x) + p2*np.sin(p1*x)+p3*x
p1, p2, p3 = 1, 0.2, 0.01
fit0=so.curve_fit(model0, xdata, ydata, p0=(p1,p2))[0]
fit1=so.curve_fit(model1, xdata, ydata, p0=(p1,p2,p3))[0]
yfit0=model0(xdata, fit0[0], fit0[1])
yfit1=model1(xdata, fit1[0], fit1[1], fit1[2])
ssq0=((yfit0-ydata)**2).sum()
ssq1=((yfit1-ydata)**2).sum()
df=len(xdata)-3
f_ratio=(ssq0-ssq1)/(ssq1/df)
p=1-ss.f.cdf(f_ratio, 1, df)
In [2]:

print f_ratio, p
0.919387419515 0.369574503394

@usethedeastar が指摘したように:剰余が正規分布している場合、非線形最小二乗 は 最大の可能性。したがって、F 検定と尤度比検定は同等です。なぜなら、 F-ratio は、尤度比 λ の単調変換です。.

または、説明的な方法については、次を参照してください。 http://www.stata.com/support/faqs/statistics/chi-squared-and-f-distributions/

他のヒント

あなたの式は私には正しく見えます。scipy.stats.norm.logpdf(x, loc=mu, scale=sigma)

と同じ結果をお勧めします。

あなたはすでにMUとSigmaの見積もりを持っているので、あなたが結果を差し込むことができる尤度比テストの機能があるとは思わない。

あなたが2つのモデルの推定値があるならば、あなたはもう一方のモデルに入れ子になっている場合、あなたは簡単に自分で計算することができます。

href="http://en.wikipedia.org/wiki/likeRihhood-ratio_test" rel="nofollow"> http://en.wikipedia.org/wiki/likeriahood-ratio_test

これは、2つの入れ子型線形モデルを比較するためのLRテストを計算するSTATSMODELSのメソッドの一部です。 https://github.com/statsmodels/statsmodels/blob./master/statsmodels/Regression/Linear_Model.py#L1531

ライセンス： CC-BY-SA と 帰属
所属していません StackOverflow
scroll top