Đừng để dữ liệu nhiễu bẻ cong đường hồi quy của bạn

Trong bài trước của loạt bài về ước lượng bền vững, chúng ta đã so sánh hồi quy OLS với các phương pháp như Huber, RANSAC, GNC-GM, GNC-TLS và ASOR để xem điều gì xảy ra khi dữ liệu nhiễm bẩn kéo lệch một đường thẳng khỏi quan hệ thực. Lần này, bài toán vẫn giữ nguyên nhưng mô hình được nới lỏng: thay vì một đường thẳng, ta cho phép đường cong.

Minh họa đường hồi quy bị kéo lệch bởi dữ liệu nhiễu

Khi đường thẳng không còn đủ

Hồi quy tuyến tính phổ biến vì đơn giản, nhanh, dễ diễn giải và thường hiệu quả bất ngờ. Nhưng quan hệ thực tế không phải lúc nào cũng chịu hợp tác với đường thẳng. Có quan hệ cong dần, có quan hệ bão hòa, có quan hệ thay đổi hành vi theo từng vùng. Khoảng cách phanh xe, phản ứng năng suất cây trồng, công suất tuabin gió hay Doppler vệ tinh đều là những ví dụ quen thuộc mà độ cong chính là một phần của bản chất quan hệ, chứ không phải thứ gây nhiễu cần bỏ qua.

Bài toán giờ đây tách thành hai câu hỏi khác nhau:

• Bộ ước lượng nên tin tưởng những quan sát nào?
• Cơ sở hàm nên có khả năng biểu diễn những hình dạng nào?

Nói cách khác: bộ ước lượng kiểm soát mức độ tin cậy, còn cơ sở hàm kiểm soát hình dạng.

Hồi quy đa thức dưới góc nhìn mô hình tuyến tính

Với N phép đo, mô hình bậc hai có dạng:

y_i = β₀ + β₁x_i + β₂x_i² + ε_i

Trong đó các hệ số β_k là tham số mô hình, còn sai số ε_i tuân theo phân phối chuẩn N(0, σ²_in). Điểm mấu chốt: y_i biến thiên phi tuyến theo x_i, nhưng mô hình vẫn tuyến tính theo các hệ số chưa biết. Hồi quy đa thức thực chất là mô hình tuyến tính trong một không gian đặc trưng đã biến đổi.

Dưới dạng ma trận: y = Φβ + ε, với Φ là ma trận thiết kế chứa các lũy thừa của x. OLS giải bài toán cực tiểu hóa chuẩn bình phương của phần dư.

Trong code, chuyển từ đường thẳng sang đa thức chỉ là thay đổi ma trận thiết kế:

def polynomial_design(x, degree):
    x = np.asarray(x, dtype=float)
    return np.vander(x, N=degree + 1, increasing=True)

Phi = polynomial_design(x, degree=3)
beta = np.linalg.lstsq(Phi, y, rcond=None)[0]

Để ổn định số học — đặc biệt khi bậc tăng cao hoặc dải giá trị x rộng — nên chuẩn hóa tâm và thang đo của biến dự báo trước khi tạo các lũy thừa. Việc chuẩn hóa giúp các số hạng bậc cao không phình quá lớn và cải thiện độ ổn định tính toán.

Các hàm cơ sở đa thức tăng dần theo bậc

Nhưng linh hoạt hơn luôn đi kèm đánh đổi. Đa thức bậc p có p+1 hệ số, nên các phương pháp như RANSAC cần mẫu tối thiểu lớn hơn khi bậc tăng. Ngoại suy cũng nhạy cảm hơn, và vì hàm cơ sở đa thức mang tính toàn cục, chỉ cần thay đổi một hệ số là có thể định hình lại toàn bộ đường cong.

Khi một đa thức toàn cục là không đủ

Spline tiếp cận theo cách khác: thay vì dùng một đa thức toàn cục, nó ghép các mảnh đa thức tại những vị trí gọi là knot (nút), đồng thời đảm bảo các mảnh liền kề nối với nhau một cách trơn tru. B-spline cung cấp một cơ sở thuận tiện cho việc này.

Mô hình spline có dạng: f(x) = Σ β_j · B_j,p(x), với B_j,p(x) là hàm cơ sở B-spline thứ j bậc p.

Cơ sở B-spline được xác định bởi bậc spline và vị trí các knot. Với một khoảng quan sát [a, b], bậc p và K knot bên trong, mỗi biên được lặp lại p+1 lần, tạo ra vector knot mở. Số hàm cơ sở thu được là M = K + p + 1. Ví dụ, spline bậc hai với 5 knot bên trong sẽ cho 8 hàm cơ sở.

Các hàm cơ sở này được xây dựng bằng đệ quy Cox–de Boor: bắt đầu từ các hàm bậc không (hằng số từng khúc), rồi nâng dần lên bậc một (tuyến tính từng khúc) và cuối cùng là bậc hai (bình phương từng khúc).

Điểm khác biệt cốt lõi so với lũy thừa của đa thức toàn cục: mỗi hàm cơ sở B-spline bậc hai chỉ có giá đỡ cục bộ — nó chỉ khác 0 trong một khoảng hữu hạn. Nhờ đó, mỗi hệ số chỉ ảnh hưởng đến một vùng giới hạn của đường cong. Với các knot bên trong đơn giản, các mảnh bình phương liền kề còn nối với nhau bằng đạo hàm bậc một liên tục.

Sau khi có các hàm cơ sở, mỗi quan sát được biểu diễn thành một hàng của ma trận thiết kế, và ta quay lại đúng mô hình quen thuộc y = Φβ + ε. Hồi quy đa thức và hồi quy spline chỉ khác nhau ở cách xây dựng Φ, còn cả hai vẫn tuyến tính theo các hệ số chưa biết khi cơ sở đã cố định.

Quy trình xây dựng hàm cơ sở B-spline bậc hai

Vấn đề ngoại lai và sáu bộ ước lượng quen thuộc

Dữ liệu thiếu có thể được xử lý đơn giản bằng cách lọc bỏ các giá trị không hữu hạn trước khi huấn luyện. Nhưng vấn đề khó hơn là những quan sát có mặt nhưng không đáng tin. Ngoại lai xuất hiện do cảm biến lỗi, nhiễu can thiệp hay điều kiện vận hành bất thường; chúng vẫn hữu hạn và thoạt nhìn hoàn toàn hợp lệ.

Điều này tạo ra bài toán "con gà – quả trứng": ta cần một mô hình đáng tin để xác định quan sát nào không nhất quán, nhưng chính những quan sát đó lại có thể bóp méo mô hình dùng để nhận diện chúng.

Tin tốt là thay đổi cơ sở hàm không buộc phải thiết kế lại bộ ước lượng bền vững. Với các bộ ước lượng dựa trên trọng số, bước hồi quy vẫn là:

β̂ = argmin Σ w_i · (y_i − φ(x_i)ᵀβ)²

Điều thay đổi chỉ là vector đặc trưng φ(x_i): với đa thức bậc ba là [1, x, x², x³], còn với spline 8 hàm là [B₁(x), …, B₈(x)].

Kết quả kiểm thử áp lực

Nhóm tác giả lặp lại cùng họ thí nghiệm bền vững cho ba mô hình đường cong — đa thức bậc hai, đa thức bậc ba và B-spline bậc hai — với sáu bộ ước lượng. Benchmark bao gồm nhiều dạng nhiễm bẩn: ngoại lai thay thế Gauss trung bình không, Gauss lệch, đều lệch, và cả quan sát sinh ra từ một quan hệ cạnh tranh mạch lạc. Tỷ lệ nhiễm bẩn quét từ 0% đến 90%, mỗi điều kiện lặp lại 30 lần Monte Carlo.

Kết quả nổi bật:

• Ở mức nhiễm bẩn Gauss lệch 30%, GNC-TLS cho NRMSE trung vị lần lượt 0,246 / 0,276 / 0,485 cho đa thức bậc hai, bậc ba và spline. ASOR có xu hướng tương tự: 0,257 / 0,307 / 0,469. Việc tăng từ 3 lên 4 hệ số chỉ gây mất mát khiêm tốn, nhưng spline 8 hàm đòi hỏi cao hơn rõ rệt.
• RANSAC phản ứng mạnh nhất. Cùng dữ liệu và mức nhiễm bẩn, NRMSE trung vị của nó tăng từ 0,273 (bậc hai) lên 0,383 (bậc ba) và 1,384 (spline). Ngay cả khi không có nhiễu, con số này đã tăng từ 0,179 lên 0,625.
• Lý do: mẫu tối thiểu của RANSAC tăng theo số hệ số (s = 3, 4, 8). Ở mức nhiễm bẩn 50% và độ tin cậy 0,999, số giả thuyết lý thuyết cần thiết tăng từ khoảng 52 lên 108 rồi 1765. Thời gian chạy trung vị cũng tăng từ khoảng 2,1 ms lên 3,4 ms và 17,4 ms.
• Mô hình giàu tham số cần nhiều dữ liệu hơn. Với spline dưới nhiễu Gauss lệch 50%, NRMSE của GNC-TLS giảm từ 1,641 ở N = 50 xuống 0,275 ở N = 5000. Trong khi đó OLS vẫn sai lệch cao dù dữ liệu tăng, vì thêm dữ liệu làm giảm phương sai chứ không loại bỏ được lực kéo hệ thống của nhiễu.
• Không có bảng xếp hạng phổ quát. Dưới nhiễu lệch, GNC-TLS và ASOR thường cho sai số dự báo thấp nhất. Nhưng khi sai số thay thế có trung bình bằng 0, khoảng cách giữa các phương pháp thu hẹp đáng kể.
• Giá đỡ cục bộ chỉ khoanh vùng méo mó, không xóa bỏ nó. Trong thí nghiệm minh họa, khi chỉ nhiễm bẩn một vùng hẹp, spline bị méo mạnh hơn bên trong vùng nhiễm (RMSE 40,9 so với 22,6 của đa thức bậc ba) nhưng nhẹ hơn ở bên ngoài (12,1 so với 13,6). Khi thêm trọng số bền vững, sai số giảm mạnh ở cả hai biểu diễn.

Giá đỡ cục bộ kiểm soát nơi tính linh hoạt của mô hình tác động. Ước lượng bền vững kiểm soát mức độ tác động của những quan sát không đáng tin. Hai cơ chế giải quyết hai bài toán khác nhau.

Một trường hợp đáng lưu ý: khi ngoại lai tạo thành một đường cong cạnh tranh mạch lạc, các bộ ước lượng bền vững vẫn bám được đường cong gốc ở mức nhiễm bẩn thấp và vừa. Nhưng khi hai cấu trúc gần cân bằng, bài toán trở nên nhập nhằng. Một bộ ước lượng đơn mô hình không có cách nào biết quan hệ nào trong hai quan hệ nhất quán nội tại mới là đúng về mặt vật lý nếu không có thêm thông tin.

Nghiên cứu tình huống 1: Đường cong công suất tuabin gió

Độ bền vững không đồng nghĩa với tính linh hoạt của mô hình. Dữ liệu SCADA của tuabin gió thường chứa các quan sát công suất thấp bất thường do cắt giảm công suất, sự cố, dừng máy hay lỗi đo lường. Những quan sát này có thể bóp méo đường cong công suất danh định.

Nhóm tác giả tách bạch hai vấn đề bằng cách trước tiên tạo một vùng mà đa thức là mô hình đúng. Dưới tốc độ định mức 12 m/s, đặt t = (v − 3)/9 và P(v) = 2000(3t² − 2t³). Vì t tuyến tính theo v, biểu thức này chính xác là một đa thức bậc ba theo tốc độ gió. Khi đó họ mô hình được chỉ định đúng, và thí nghiệm chỉ còn đo tác động của nhiễu: RMSE dự báo khoảng 91,4 kW với OLS và chỉ 2,7 kW với ASOR.

Sau đó mở rộng miền quan sát vượt tốc độ định mức, thêm vùng công suất bão hòa 2000 kW từ 12 đến 14 m/s. Lúc này quan hệ trở thành từng khúc: vùng tăng dưới định mức nối với vùng bão hòa. Không đa thức hữu hạn bậc khác hằng nào có thể tái tạo chính xác hình dạng này — nếu một đa thức bằng 2000 trên toàn khoảng [12, 14] thì nó phải bằng 2000 ở mọi nơi.

Tăng bậc đa thức lên bậc bốn giúp cải thiện, nhưng cả năm hàm cơ sở vẫn mang tính toàn cục: cùng một bộ hệ số phải đồng thời mô tả vùng tăng chậm, vùng dốc giữa, chuyển tiếp quanh tốc độ định mức và vùng bão hòa. Trên cùng một lần minh họa, đa thức bậc bốn với ASOR cho RMSE khoảng 14,5 kW, so với 4,2 kW của B-spline bậc hai có năm knot bên trong.

Ước lượng bền vững có thể ngăn quan sát bất thường bẻ cong một mô hình tốt. Nó không thể biến một cơ sở hàm sai thành đúng.

Nghiên cứu tình huống 2: Dự đoán Doppler GNSS

Doppler GNSS biến thiên trơn tru theo chuyển động tương đối giữa vệ tinh và máy thu. Dự đoán ngắn hạn giúp lấp khoảng trống đo đạc ngắn và cung cấp xu hướng Doppler kỳ vọng để so sánh với các phép đo bất thường do nhiễu hay mất bám tín hiệu.

Quy trình là: lấy mẫu Doppler gần nhất → khớp đa thức cục bộ → dự đoán ngắn hạn. Trong ví dụ mô phỏng có cả mẫu thiếu lẫn mẫu bị hỏng, OLS cố giải thích các phép đo hỏng nên đường khớp bị bẻ về phía cụm dữ liệu lỗi, trong khi ASOR giảm ảnh hưởng của chúng.

Điểm quan trọng là nội suy và ngoại suy không phải cùng một bài toán. Một sai số vừa phải trong vùng quan sát có thể phình to đáng kể khi ngoại suy. Trong minh họa này, RMSE dự báo ngắn hạn khoảng 46,4 Hz với OLS và 7,9 Hz với ASOR. Tuy nhiên, tác giả nhấn mạnh lập luận mang tính cục bộ: không giả định quỹ đạo Doppler tuân theo một đa thức bậc ba vô hạn, mà chỉ dùng đa thức bậc thấp như xấp xỉ gọn trong một cửa sổ quan sát ngắn.

Điều gì còn lại khi đường thẳng hóa thành đường cong?

Bài học hữu ích nhất không gắn với một bậc đa thức hay một loại spline cụ thể, mà nằm ở góc nhìn cơ sở hàm. Đường thẳng dùng [1, x], đa thức dùng [1, x, x², …, x^p], spline thay lũy thừa toàn cục bằng các hàm có giá đỡ cục bộ. Khi cơ sở đã cố định, tất cả đều quay về dạng y = Φβ + ε. Đó là lý do các bộ ước lượng bền vững chuyển sang tự nhiên đến vậy.

Quy trình thực hành rút ra:

• Chọn cơ sở hàm có khả năng biểu diễn quan hệ cần quan tâm.
• Chuẩn hóa tâm và thang đo biến dự báo khi lũy thừa đa thức trở nên kém ổn định.
• Loại bỏ các hàng thiếu có thể nhận diện trước khi khớp, nếu phân tích theo trường hợp đầy đủ là phù hợp.
• Dùng ước lượng bền vững khi các quan sát hữu hạn vẫn có thể bị nhiễm bẩn.
• Dùng đủ số quan sát để hỗ trợ cơ sở hàm đã chọn.
• Thận trọng hơn với ngoại suy so với nội suy.

Minh họa quy trình chọn cơ sở hàm và bộ ước lượng bền vững

Với độc giả làm việc trong lĩnh vực khoa học dữ liệu tại Việt Nam — từ phân tích dữ liệu cảm biến IoT trong nông nghiệp, giám sát thiết bị công nghiệp đến các hệ thống định vị — bài học này khá thực dụng: đừng kỳ vọng một thuật toán "chống nhiễu" có thể cứu vãn một lựa chọn mô hình sai. Cơ sở hàm giàu biểu đạt hơn giúp mô hình diễn đạt tốt hơn, nhưng không làm cho quan sát không đáng tin trở nên đáng tin. Ước lượng bền vững giúp mô hình ít nhạy cảm với nhiễu, nhưng không cứu được một mô hình không phù hợp.

Đó là hai bài toán khác nhau — và nhầm lẫn giữa chúng là nguồn gốc của không ít đường cong bị bẻ cong một cách khó hiểu trong thực tế.