置信區間與最小二乘擬合python

Question

如何計算最小二乘擬合（scipy.optimize.leastsq）在Python中的置信區間？

Answer 1

我會使用bootstrapping方法。
在這裏看到：在嘈雜的高斯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

我不知道你所說的置信區間的意思。

一般而言，leastsq對您嘗試最小化的函數並不瞭解太多，所以它不能真正給出置信區間。但是，它確實返回了Hessian的估計值，換句話說就是將第二個導數推廣到多維問題。

正如函數的文檔字符串中所暗示的那樣，您可以使用該信息以及殘差（擬合解與實際數據之間的差異）來計算參數估計的協方差，這是對置信度的局部猜測間隔。

請注意，它只是一個本地信息，我懷疑，只有當您的目標函數嚴格凸時，您纔可以嚴格地得出結論。我沒有任何證據或對該聲明的引用:)。

Answer 3

估計置信區間（CI）的最簡單方法是將標準誤差（標準偏差）乘以常數。要計算常數，您需要知道自由度（DOF）的數量和要計算CI的置信度。以這種方式估計的CI有時稱爲漸近CI。 您可以在Motulsky & Christopoulos（google books）的「使用線性和非線性迴歸擬合模型生物數據」中詳細瞭解它。同一本書（或非常相似）免費提供as a manual for author's software。您也可以閱讀how to calculate CI using the C++ Boost.Math library。在這個例子中，CI是針對一個變量的分佈進行計算的。在最小二乘擬合的情況下，DOF不是N -1，而是N-M，其中M是參數的數量。在Python中執行相同的操作應該很容易。

這是最簡單的估計。我不知道zephyr提出的bootstrapping方法，但它可能比我寫的方法更可靠。