2011-08-19 32 views
7

Tôi đã thực hiện một số mô phỏng vật lý Monte Carlo với Python và tôi không thể xác định sai số chuẩn cho các hệ số của một hình vuông không tuyến tính nhỏ nhất.Lỗi chuẩn trong hồi quy phi tuyến

Ban đầu, tôi đã sử dụng SciPy scipy.stats.linregress cho mô hình của mình vì tôi nghĩ nó sẽ là mô hình tuyến tính nhưng nhận thấy nó thực sự là một loại chức năng điện. Sau đó tôi đã sử dụng số polyfit của NumPy với mức độ tự do là 2 nhưng tôi không thể tìm thấy cách nào để xác định lỗi chuẩn của các hệ số.

Tôi biết gnuplot có thể xác định lỗi cho tôi nhưng tôi cần phải làm phù hợp với hơn 30 trường hợp khác nhau. Tôi đã tự hỏi nếu có ai biết anyway cho Python để đọc các lỗi tiêu chuẩn từ gnuplot hoặc là có một số thư viện khác tôi có thể sử dụng?

Trả lời

2

nó trông giống như gnuplot sử dụng levenberg-marquardt và có một con trăn implementation available - bạn có thể nhận được ước tính lỗi từ các thuộc tính mpfit.covar (tình cờ, bạn nên lo lắng về những gì dự toán lỗi "có nghĩa là" - được các thông số khác được phép điều chỉnh để bù đắp, ví dụ?)

+0

Cảm ơn bạn đã liên kết! Cuối cùng tôi đã không sử dụng mpfit nhưng tài liệu đã giúp tôi rất nhiều để hiểu curve_fit cho scipy! – syntaxing

10

Cuối cùng tìm thấy câu trả lời cho câu hỏi dài này! Tôi hy vọng điều này có thể ít nhất là tiết kiệm cho ai đó một vài giờ nghiên cứu vô vọng cho chủ đề này. Scipy có một chức năng đặc biệt gọi là curve_fit trong phần tối ưu hóa của nó. Nó sử dụng phương pháp tối thiểu để xác định các hệ số và tốt nhất của tất cả, nó mang lại cho bạn ma trận hiệp phương sai. Ma trận hiệp phương sai chứa phương sai của mỗi hệ số. Chính xác hơn, đường chéo của ma trận là phương sai và bằng cách căn lề các giá trị, lỗi chuẩn của mỗi hệ số có thể được xác định! Scipy không có nhiều tài liệu hướng dẫn cho điều này vì vậy đây là một mẫu mã cho một sự hiểu biết tốt hơn:

import numpy as np 
from scipy.optimize import curve_fit 
import matplotlib.pyplot as plot 


def func(x,a,b,c): 
    return a*x**2 + b*x + C#Refer [1] 

x = np.linspace(0,4,50) 
y = func(x,2.6,2,3) + 4*np.random.normal(size=len(x)) #Refer [2] 


coeff, var_matrix = curve_fit(func,x,y) 
variance = np.diagonal(var_matrix) #Refer [3] 

SE = np.sqrt(variance) #Refer [4] 

#======Making a dictionary to print results======== 
results = {'a':[coeff[0],SE[0]],'b':[coeff[1],SE[1]],'c':[coeff[2],SE[2]]} 

print "Coeff\tValue\t\tError" 
for v,c in results.iteritems(): 
    print v,"\t",c[0],"\t",c[1] 
#========End Results Printing================= 

y2 = func(x,coeff[0],coeff[1],coeff[2]) #Saves the y values for the fitted model 

plot.plot(x,y) 
plot.plot(x,y2) 

plot.show() 
  1. gì hàm này trả về là rất quan trọng vì nó xác định những gì sẽ được sử dụng để phù hợp cho mô hình
  2. Sử dụng chức năng để tạo ra một số dữ liệu tùy ý + một số tiếng ồn
  3. Tiết kiệm ma trận hiệp phương sai của đường chéo để một ma trận 1D mà chỉ là một mảng bình thường
  4. Quảng trường rễ phương sai để có được những sai số chuẩn (SE)
Các vấn đề liên quan