Hướng dẫn weighted nonlinear regression python - python hồi quy phi tuyến tính có trọng số
|
Một trong những ứng dụng chính của bình phương tối thiểu phi tuyến là hồi quy phi tuyến hoặc phù hợp với đường cong. Đó là bởi các cặp đã cho $ \ left \ {(t_i, y_i) \: i = 1, \ ldots, n \ right \} $ ước tính tham số $ \ mathbf {x} $ xác định hàm phi tuyến $ \ varphi (t; Mathbf {x}) $, giả sử mô hình: \ started {phương trình} y_i = \ varphi (t_i; Trong đó $ \ epsilon_i $ là lỗi đo lường (quan sát). Trong ước tính bình phương nhỏ nhất, chúng tôi tìm kiếm $ \ mathbf {x} $ là giải pháp của vấn đề tối ưu hóa sau: \ start T_I; Công thức như vậy là trực quan và thuyết phục theo quan điểm toán học. Từ quan điểm xác suất, giải pháp bình phương nhỏ nhất được biết đến là ước tính khả năng tối đa, với điều kiện là tất cả $ \ epsilon_i $ đều độc lập và thường được phân phối các biến ngẫu nhiên. Vì vậy, về mặt lý thuyết, nó không tối ưu khi $ \ epsilon_i $ có phân phối khác với bình thường. Mặc dù trong thực tiễn kỹ thuật, nó không quan trọng, tức là nếu các lỗi hoạt động như một số biến ngẫu nhiên hợp lý với 0 có nghĩa là kết quả của ước tính bình phương nhỏ nhất sẽ là thỏa đáng. Các vấn đề thực sự bắt đầu khi dữ liệu bị ô nhiễm bởi các ngoại lệ (các phép đo hoàn toàn sai). Trong trường hợp này, giải pháp bình phương nhỏ nhất có thể trở nên thiên vị đáng kể để tránh dư rất cao trên các ngoại lệ. Để giải thích một cách định tính lý do tại sao nó đang xảy ra, hãy xem xét vấn đề bình phương nhỏ nhất đơn giản nhất: \ Bắt đầu {Align} \ frac {1} {2} \ sum_ {i = 1}^n ( \ end {align} và giải pháp là giá trị trung bình: \ start Bây giờ hãy tưởng tượng: $ a = 1, n = 4, x_1 = 0.9, x_2 = 1.05, x_3 = 0,95, x_4 = 10 \ text {(ngoại lệ)} $. Giải pháp $ a^* = 3.225 $, tức là bị hủy hoại hoàn toàn bởi một ngoại lệ duy nhất. Một trong những công cụ ước tính mạnh mẽ nổi tiếng là người ước tính L1, trong đó tổng các giá trị tuyệt đối của phần dư được giảm thiểu. Để trình diễn, một lần nữa, hãy xem xét vấn đề đơn giản nhất: \ Bắt đầu {Công thức} \ sum_ {i = 1}^n | a - x_i | \ Rightarrow \ min_a \ end {phương trình} Và giải pháp nổi tiếng là trung bình của $ \ {x_i \} $, sao cho số lượng nhỏ các ngoại lệ sẽ không ảnh hưởng đến giải pháp. Trong vấn đề đồ chơi của chúng tôi, chúng tôi có $ a^* = 1 $. Nhược điểm duy nhất của người ước tính L1 là vấn đề tối ưu hóa phát sinh là khó khăn, vì chức năng không phân biệt ở mọi nơi, điều này đặc biệt rắc rối để tối ưu hóa phi tuyến hiệu quả. Điều đó có nghĩa là chúng ta tốt hơn để ở lại với các vấn đề khác nhau, nhưng bằng cách nào đó kết hợp sự mạnh mẽ trong ước tính. Để thực hiện điều này, chúng tôi giới thiệu một hàm sublinear $ \ rho (z) $ (nghĩa là sự tăng trưởng của nó phải chậm hơn tuyến tính) và xây dựng một vấn đề tối ưu hóa giống như bình phương mới nhất: \ Begin {phương trình} \ frac {1} {2} \ sum_ {i = 1}^n \ rho \ left ((\ varphi (t_i; x) - y_i)^2 \ right) Hóa ra vấn đề này có thể được giảm xuống thành bình phương tối thiểu phi tuyến tiêu chuẩn bằng cách sửa đổi một vectơ của phần dư và ma trận Jacobian trên mỗi lần lặp, do đó tính toán độ dốc và xấp xỉ Hessian phù hợp với hàm của hàm mục tiêu. Tham khảo bài báo để biết chi tiết. Các lựa chọn được triển khai của $ \ rho $ là như sau: Các hàm 2 và 3 tương đối nhẹ và gây mất giá trị xấp xỉ tuyệt đối cho phần dư lớn. Hai chức năng cuối cùng là thăng hoa mạnh mẽ và suy giảm đáng kể cho các ngoại lệ. Các hàm mất ở trên được viết với giả định rằng ngưỡng mềm giữa các inlin và ngoại lệ bằng 1,0. Để khái quát hóa, chúng tôi giới thiệu tham số tỷ lệ $ c $ và đánh giá tổn thất là \ Begin {phương trình} \ hat {\ rho} (r^2) = c^2 \ rho \ trái (\ trái (\ frac {r} { C} \ right)^2 \ right) \ end {phương trình} Lưu ý rằng nếu phần dư nhỏ, chúng ta có $ \ hat {\ rho} (r^2) \ xấp xỉ \ rho (r^2) \ xấp xỉ r^2 $ cho bất kỳ $ \ rho (z) $ được xác định ở trên. Để minh họa hành vi của các chức năng được giới thiệu, chúng tôi xây dựng một lô: In [1]: import numpy as np %matplotlib inline import matplotlib.pyplot as plt from matplotlib import rcParams rcParams['figure.figsize'] = (10, 6) rcParams['legend.fontsize'] = 16 rcParams['axes.labelsize'] = 16 In [2]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2) In [3]: plt.plot(r, linear, label='linear') plt.plot(r, huber, label='huber') plt.plot(r, soft_l1, label='soft_l1') plt.plot(r, cauchy, label='cauchy') plt.plot(r, arctan, label='arctan') plt.xlabel("$r$") plt.ylabel(r"$\rho(r^2)$") plt.legend(loc='upper left'); Out[3]: Bây giờ chúng tôi sẽ chỉ ra cách các chức năng tổn thất mạnh mẽ hoạt động trên một ví dụ mô hình. Chúng tôi xác định hàm mô hình là \ start Có thể mô hình hóa sự dịch chuyển quan sát của bộ tạo dao động ẩm tuyến tính. In [4]: def generate_data(t, A, sigma, omega, noise=0, n_outliers=0, random_state=0): y = A * np.exp(-sigma * t) * np.sin(omega * t) rnd = np.random.RandomState(random_state) error = noise * rnd.randn(t.size) outliers = rnd.randint(0, t.size, n_outliers) error[outliers] *= 35 return y + error In [5]: A = 2 sigma = 0.1 omega = 0.1 * 2 * np.pi x_true = np.array([A, sigma, omega]) noise = 0.1 t_min = 0 t_max = 30 Dữ liệu để lắp các tham số sẽ chứa 3 ngoại lệ: In [6]: t_train = np.linspace(t_min, t_max, 30) y_train = generate_data(t_train, A, sigma, omega, noise=noise, n_outliers=4) Xác định chức năng Phần dư tính toán để giảm thiểu bình phương nhỏ nhất: In [7]: def fun(x, t, y): return x[0] * np.exp(-x[1] * t) * np.sin(x[2] * t) - y Sử dụng tất cả những cái như ước tính ban đầu. In [9]: from scipy.optimize import least_squares Chạy bình phương tối thiểu tiêu chuẩn: In [10]: res_lsq = least_squares(fun, x0, args=(t_train, y_train)) Chạy bình phương tối thiểu mạnh mẽ với r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)5, đặt r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)6 thành 0,1, có nghĩa là phần dư inliter thấp hơn xấp xỉ 0,1. In [11]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)0 Xác định dữ liệu để vẽ các đường cong đầy đủ. In [12]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)1 Dự đoán tính toán với các tham số tìm thấy: In [13]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)2 In [14]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)3 Out[14]: r = np.linspace(0, 5, 100) linear = r**2 huber = r**2 huber[huber > 1] = 2 * r[huber > 1] - 1 soft_l1 = 2 * (np.sqrt(1 + r**2) - 1) cauchy = np.log1p(r**2) arctan = np.arctan(r**2)4 Chúng tôi thấy rõ rằng bình phương tối thiểu tiêu chuẩn phản ứng đáng kể trên các ngoại lệ và đưa ra giải pháp rất thiên vị, trong khi các bình phương tối thiểu mạnh mẽ (với L1-Loss mềm) cung cấp giải pháp rất gần với mô hình thực. Phần tác giả: Thị trưởng Nikolay |
