python - scipy pythonでのleastsqフィットによる信頼区間

Question

Pythonで最小二乗適合(scipy.optimize.leastsq)の信頼区間を計算するには?

Answer 1

ブートストラップ法を使用します。
ここを参照してください: http://phe.rockefeller.edu/LogletLab/whitepaper/node17.html

ノイズの多いガウス分布の簡単な例:

x = arange(-10, 10, 0.01)

# model function
def f(p):
    mu, s = p
    return exp(-(x-mu)**2/(2*s**2))

# create error function for dataset    
def fff(d):
    def ff(p):
        return d-f(p)
    return ff

# create noisy dataset from model
def noisy_data(p):
    return f(p)+normal(0,0.1,len(x))

# fit dataset to model with least squares    
def fit(d):
    ff = fff(d)
    p = leastsq(ff,[0,1])[0]
    return p

# bootstrap estimation        
def bootstrap(d):
    p0 = fit(d)
    residuals = f(p0)-d
    s_residuals = std(residuals)

    ps = []
    for i in range(1000):
        new_d = d+normal(0,s_residuals,len(d))
        ps.append(fit(new_d))

    ps = array(ps)
    mean_params = mean(ps,0)
    std_params = std(ps,0)

    return mean_params, std_params

data = noisy_data([0.5, 2.1])
mean_params, std_params = bootstrap(data)

print "95% confidence interval:"
print "mu: ", mean_params[0], " +/- ", std_params[0]*1.95996
print "sigma: ", mean_params[1], " +/- ", std_params[1]*1.95996

Answer 2

信頼区間 (CI) を推定する最も簡単な方法は、標準誤差 (標準偏差) に定数を掛けることです。定数を計算するには、自由度 (DOF) の数と、CI を計算する信頼レベルを知る必要があります。この方法で推定された CI は、漸近的 CI と呼ばれることがあります。詳細については、Motulsky & Christopoulos ( Google ブックス) による「線形および非線形回帰を使用した生物学的データへのモデルの適合」を参照してください。同じ本 (または非常に類似した本) が、著者のソフトウェアのマニュアルとして無料で入手できます。

C++ Boost.Math ライブラリを使用して CI を計算する方法もお読みください。この例では、CI は 1 つの変数の分布に対して計算されます。最小二乗フィッティングの場合、DOF はN -1 ではなくNMです。ここで、Mはパラメーターの数です。Python で同じことを行うのは簡単なはずです。

これが最も単純な見積もりです。zephyrが提案しているブートストラップ方法は知りませんが、私が書いた方法よりは信頼できるかもしれません。