2014-11-25 16 views
6

Tôi có một câu hỏi mà tôi đã đấu tranh trong nhiều ngày với bây giờ.Tính toán dải tin cậy của hình vuông vừa vặn nhỏ nhất

Làm cách nào để tính toán dải tin cậy (95%) phù hợp?

đường cong Lắp dữ liệu là công việc hàng ngày của tất cả các nhà vật lý - vì vậy tôi nghĩ rằng điều này nên được thực hiện ở đâu đó - nhưng tôi không thể tìm thấy một thực hiện cho điều này không làm tôi biết làm thế nào để làm điều này về mặt toán học .

Điều duy nhất tôi tìm thấy là seaborn thực hiện công việc tốt đẹp cho tuyến tính nhỏ nhất vuông.

import numpy as np                                                       
from matplotlib import pyplot as plt 
import seaborn as sns 
import pandas as pd 

x = np.linspace(0,10) 
y = 3*np.random.randn(50) + x 

data = {'x':x, 'y':y} 
frame = pd.DataFrame(data, columns=['x', 'y']) 
sns.lmplot('x', 'y', frame, ci=95) 

plt.savefig("confidence_band.pdf") 

linear-least-square

Nhưng điều này chỉ là tuyến tính ít nhất vuông. Khi tôi muốn phù hợp với ví dụ: một đường cong bão hòa như saturation-eqn, tôi hơi say.

Chắc chắn, tôi có thể tính toán phân phối t từ lỗi std của phương pháp tối thiểu như scipy.optimize.curve_fit nhưng đó không phải là những gì tôi đang tìm kiếm.

Cảm ơn bạn đã trợ giúp !!

Trả lời

6

Bạn có thể thực hiện điều này dễ dàng bằng mô-đun StatsModels.

Đồng thời xem this examplethis answer.

Dưới đây là một câu trả lời cho câu hỏi của bạn:

import numpy as np 
from matplotlib import pyplot as plt 

import statsmodels.api as sm 
from statsmodels.stats.outliers_influence import summary_table 

x = np.linspace(0,10) 
y = 3*np.random.randn(50) + x 
X = sm.add_constant(x) 
res = sm.OLS(y, X).fit() 

st, data, ss2 = summary_table(res, alpha=0.05) 
fittedvalues = data[:,2] 
predict_mean_se = data[:,3] 
predict_mean_ci_low, predict_mean_ci_upp = data[:,4:6].T 
predict_ci_low, predict_ci_upp = data[:,6:8].T 

fig, ax = plt.subplots(figsize=(8,6)) 
ax.plot(x, y, 'o', label="data") 
ax.plot(X, fittedvalues, 'r-', label='OLS') 
ax.plot(X, predict_ci_low, 'b--') 
ax.plot(X, predict_ci_upp, 'b--') 
ax.plot(X, predict_mean_ci_low, 'g--') 
ax.plot(X, predict_mean_ci_upp, 'g--') 
ax.legend(loc='best'); 
plt.show() 

Example

+0

Thật không may, điều này hiện chỉ có sẵn trong statsmodels cho các chức năng tuyến tính, và sẽ có sẵn cho các mô hình tuyến tính tổng quát hóa trong phiên bản tiếp theo, nhưng chưa cho các chức năng phi tuyến tính tổng quát. – user333700

0

kmpfit 's confidence_band() tính toán biên độ tín nhiệm đối với phi tuyến tính bình phương nhỏ nhất. Vào đây để đường cong bão hòa của bạn:

from pylab import * 
from kapteyn import kmpfit 

def model(p, x): 
    a, b = p 
    return a*(1-np.exp(b*x)) 

x = np.linspace(0, 10, 100) 
y = .1*np.random.randn(x.size) + model([1, -.4], x) 

fit = kmpfit.simplefit(model, [.1, -.1], x, y) 
a, b = fit.params 
dfdp = [1-np.exp(b*x), -a*x*np.exp(b*x)] 
yhat, upper, lower = fit.confidence_band(x, dfdp, 0.95, model) 

scatter(x, y, marker='.', color='#0000ba') 
for i, l in enumerate((upper, lower, yhat)): 
    plot(x, l, c='g' if i == 2 else 'r', lw=2) 
savefig('kmpfit confidence bands.png', bbox_inches='tight') 

Các dfdp là các đạo hàm riêng ∂f/∂p của mô hình f = a * (1-e^(b * x)) Đối với mỗi tham số p với (ví dụ: , a và b), xem số answer của tôi cho một câu hỏi tương tự cho các liên kết nền. Và đây là kết quả:

kmpfit confidence bands

Các vấn đề liên quan