2013-01-16 55 views
6

Tôi gặp một chút rắc rối khi lắp đường cong vào một số dữ liệu, nhưng không thể tìm ra nơi tôi đang đi sai.Đường cong phân rã theo cấp số nhân phù hợp với vũng nước và scipy

Trong quá khứ tôi đã làm điều này với numpy.linalg.lstsq cho hàm mũ và scipy.optimize.curve_fit cho các chức năng sigmoid. Lần này tôi muốn tạo một kịch bản lệnh cho phép tôi chỉ định các hàm khác nhau, xác định các tham số và kiểm tra sự phù hợp của chúng với dữ liệu. Trong khi làm điều này tôi nhận thấy rằng Scipy leastsq và Numpy lstsq dường như cung cấp các câu trả lời khác nhau cho cùng một bộ dữ liệu và cùng một chức năng. Hàm này chỉ đơn giản là y = e^(l*x) và bị hạn chế sao cho y=1 tại x=0.

Đường xu hướng Excel đồng ý với kết quả Numpy lstsq, nhưng khi Scipy leastsq có thể thực hiện bất kỳ chức năng nào, bạn nên tìm ra vấn đề là gì.

import scipy.optimize as optimize 
import numpy as np 
import matplotlib.pyplot as plt 

## Sampled data 
x = np.array([0, 14, 37, 975, 2013, 2095, 2147]) 
y = np.array([1.0, 0.764317544, 0.647136491, 0.070803763, 0.003630962,  0.001485394,  0.000495131]) 

# function 
fp = lambda p, x: np.exp(p*x) 

# error function 
e = lambda p, x, y: (fp(p, x) - y) 

# using scipy least squares 
l1, s = optimize.leastsq(e, -0.004, args=(x,y)) 
print l1 
# [-0.0132281] 


# using numpy least squares 
l2 = np.linalg.lstsq(np.vstack([x, np.zeros(len(x))]).T,np.log(y))[0][0] 
print l2 
# -0.00313461628963 (same answer as Excel trend line) 

# smooth x for plotting 
x_ = np.arange(0, x[-1], 0.2) 

plt.figure() 
plt.plot(x, y, 'rx', x_, fp(l1, x_), 'b-', x_, fp(l2, x_), 'g-') 
plt.show() 

Edit - thêm thông tin

Các MWe trên đã bao gồm một mẫu nhỏ của tập dữ liệu. Khi lắp dữ liệu thực tế, đường cong scipy.optimize.curve_fit trình bày R^2 là 0,82, trong khi đường cong numpy.linalg.lstsq, tương tự như được tính toán bằng Excel, có R^2 trong tổng số 0,41 .

Trả lời

4

Bạn đang giảm thiểu các chức năng lỗi khác nhau.

Khi bạn sử dụng numpy.linalg.lstsq, hàm lỗi được giảm thiểu là

np.sum((np.log(y) - p * x)**2) 

khi scipy.optimize.leastsq giảm thiểu chức năng

np.sum((y - np.exp(p * x))**2) 

Các trường hợp đầu tiên đòi hỏi một sự phụ thuộc tuyến tính giữa các biến phụ thuộc và độc lập, nhưng giải pháp được biết đến về mặt hậu môn, trong khi giải pháp thứ hai có thể xử lý bất kỳ sự phụ thuộc nào, nhưng dựa trên một phương pháp lặp lại.

Trên một lưu ý riêng, tôi không thể kiểm tra nó ngay bây giờ, nhưng khi sử dụng numpy.linalg.lstsq, tôi bạn không cần phải vstack một dãy số không, các công trình sau đây cũng như:

l2 = np.linalg.lstsq(x[:, None], np.log(y))[0][0] 
+0

Cảm ơn @ Jaime - câu trả lời tuyệt vời!Thật không may là kiến ​​thức toán học của tôi không phải là tuyệt vời; là viết hay sai [cũng thấy sửa đổi ở trên], hoặc chúng chỉ khác về cơ bản ...? Các hàm ý cho các hàm khác, ví dụ, nếu tôi muốn kiểm tra sự phù hợp của đường cong Sigmoid hoặc Gompertz với cùng một dữ liệu? – StacyR

+0

@StacyR Tôi không có kiến ​​thức để trả lời đúng câu hỏi của bạn, nhưng tôi khá chắc chắn rằng phù hợp với một hàm mũ như bạn đã làm với 'np.linalg.lstsq' chỉ là một mẹo nhanh chóng mà không tính toán lỗi đúng cách. Có một số cuộc thảo luận (khó khăn cho tôi để làm theo) ở đây: http://mathworld.wolfram.com/LeastSquaresFittingExponential.html Nếu bạn không muốn lặn sâu vào công cụ này, tôi sẽ đi với phương pháp của scipy cho tất cả mọi thứ: nó nên cung cấp cho phù hợp hơn, và kết quả của bạn sẽ phù hợp cho tất cả các chức năng. – Jaime

+0

cảm ơn một lần nữa! Tôi đã làm một số nghiên cứu thêm về điều này và, như bạn đã đề cập, đã tìm thấy rằng phương pháp 'np.linalg.lstsq' quá trọng số lỗi y ở các giá trị x thấp. Liên kết mà bạn đã chia sẻ và một số tài nguyên khác tôi tìm thấy, cho phép tôi lấy được một phương pháp phân tích khác (điều làm cho nó phức tạp là ràng buộc --- tất cả các sách mô tả phương pháp cho y = a * e^b * x thay vì hơn y = e^b * x), tuy nhiên, điều này cũng tạo ra một đường cong phù hợp tồi tệ hơn so với 'scipy.optimize.leastsq' lặp lại. – StacyR

1

Để giải thích một chút về điểm Jaime, bất kỳ sự chuyển đổi phi tuyến tính nào của dữ liệu sẽ dẫn đến một hàm lỗi khác nhau và do đó dẫn đến các giải pháp khác nhau. Những điều này sẽ dẫn đến khoảng tin cậy khác nhau cho các thông số phù hợp. Vì vậy, bạn có ba tiêu chí có thể sử dụng để đưa ra quyết định: lỗi nào bạn muốn giảm thiểu, thông số nào bạn muốn tự tin hơn và cuối cùng, nếu bạn đang sử dụng phụ kiện để dự đoán giá trị, phương pháp nào mang lại ít lỗi hơn giá trị dự đoán. Chơi xung quanh một chút về mặt phân tích và trong Excel gợi ý rằng các loại nhiễu khác nhau trong dữ liệu (ví dụ: nếu chức năng tiếng ồn quy mô biên độ, ảnh hưởng đến thời gian không đổi hoặc phụ gia) dẫn đến các lựa chọn giải pháp khác nhau.

Tôi cũng sẽ thêm rằng trong khi thủ thuật này "hoạt động" cho phân rã theo hàm mũ thành 0, nó không thể được sử dụng trong trường hợp tổng quát (và phổ biến) của số mũ bị giảm (tăng hoặc giảm) thành giá trị không thể được giả định là 0.

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