Phát triển chương trình con làm khớp dữ liệu với nhiều mô hình

Nội dung text: Phát triển chương trình con làm khớp dữ liệu với nhiều mô hình

1. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 PHÁT TRIỂN CHƯƠNG TRÌNH CON LÀM KHỚP DỮ LIỆU VỚI NHIỀU MÔ HÌNH ThS. Nguyễn Ngọc Anh1 ThS. Trương Văn Minh2 TÓM TẮT Hiện nay, có rất nhiều phần mềm máy tính cho phép người dùng làm khớp dữ liệu thực nghiệm với dạng hàm tùy ý nhập bởi người dùng. Tuy nhiên, các chương trình này có dạng đóng (đối với các chương trình thương mại) hoặc có hệ thống thư viện liên kết rất phức tạp (đối với các chương trình mã nguồn mở). Do đó, việc tận dụng thư viện của các chương trình này để nhúng vào các chương trình phần mềm nhỏ tự thiết kế là không thích hợp. Bài báo này đưa ra bộ chương trình con, cho phép người dùng làm khớp số liệu thực nghiệm với dạng hàm tùy ý, được viết bằng ngôn ngữ C++, có cấu trúc đơn giản, gói gọn trong một tập tin chỉ dài 438 dòng, thuận tiện để nhúng vào các chương trình tự phát triển. Kết quả thu được bằng chương trình được so sánh với ROOT. Từ khóa:Chương trình làm khớp nền C++, thuật toán làm khớp Levenberg– Marquardt 1. Giới thiệu trình quá nặng; hệ điều hành không hỗ Làm khớp dữ liệu theo một mô trợ; Khi đó các phần mềm tự viết sẽ hình (dạng hàm) là một thủ tục được tiến là một giải pháp thích hợp. hành rất phổ biến trong phân tích số liệu Bộ chương trình con được cung (phân tích phổ, xây dựng mô hình, xác cấp trong bài báo này cho phép người định các tham số để nội suy, ngoại suy). dùng nhúng vào trong các phần mềm tự Các thủ tục này có thể được thực hiện bởi viết, để thực thi tác vụ làm khớp số liệu các chương trình có giao diện trực quan theo mô hình bất kỳ do người dùng khai như Origin [1], SciDavis [2] hoặc các báo, sử dụng thuật toán LEVENBERG- chương trình dưới dạng lệnh thực thi như MARQUARDT [7]. Chương trình cho ROOT [3], R [4], Matlab [5], Gnuplot phép người dùng lựa chọn làm khớp có [6]. Tuy nhiên, một số là các chương trọng số hoặc không có trọng số. Bộ trình thương mại (Origin, Matlab), do đó chương trình con này có kích thước rất người sử dụng sẽ phải bỏ ra một chi phí nhỏ, chỉ ~12 kb, gói gọn trong một tập không nhỏ để trang bị phần mềm. Tiếp tin *.h, thuận tiện để người dùng khai nữa, các chương trình này thường có bộ báo trong chương trình chính. Ngôn thư viện đi kèm rất lớn, và liên kết với ngữ được sử dụng là C++.Biên dịch nhau rất phức tạp. Do, đó việc nhúng các bằng GNU g++ [8]. thư viện này vào các chương trình nhỏ tự Bộ chương trình con được hiệu viết là rất phức tạp, và làm tăng kích lực hóa bằng cách so sánh kết quả với thước của chương trình. chương trình mã nguồn mở đã được Trong thực tế, tùy thuộc vào tình chứng nhận và sử dụng rộng rãi trên các huống cụ thể, việc sử dụng các phần phòng thí nghiệm trên thế giới, mềm lớn kể trên để làm khớp không ROOT.Trong báo cáo này, bộ số liệu phải lúc nào cũng thuận lợi: chương đã được sử dụng để so sánh. 1Viện Nghiên cứu Hạt nhân Đà Lạt 2 Trường Đại học Đồng Nai 122
2. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 2. Thuật toán v s tr (6 ) Thuật toán LEVENBERG- MARQUARDT Xét bộ số liệu với n điểm thực là vectơ biến thiên của vectơ nghiệm (Xi,Yi), mô hình cần làm khớp tham số . là F(X,α), với α là vectơ tham số {α1,α2,α3, .,αm). Theo đó: Giải phương trình (3) cho phép Yi=F(Xi,α1,α2,α3, .,αm) (1) xác định , từ đó xác định được mới. Thủ tục này lặp đi lặp lại nhiều lần Để xác định các tham số tự do, ta sử dụng phương pháp bình phương cho tới khi hội tụ. Phương pháp tối thiểu [9]. Phương pháp này đòi hỏi LEVENBERG-MARQUARDT, bổ phải xác định αsao cho là cực tiểu: sung thêm vào thuật toán 2 tham số và , nhằm cải thiện khả năng hội tụ của (2) quá trình khớp. Thuật toán có thể được mô tả ngắn gọn, từng bước một như sau, lưu Trong đó là trọng số tương ứng đồ thuật toán được đưa ra trong Hình 1: với điểm số liệu thứ i. cực tiểu khi: 1. Đặt , , n=0. 2. Xác định từ phương trình: (3) (7) Đối với các hàm tuyến tính, hệ m phương trình nói trên có thể được giải ra nghiệm xác định bằng phương Với , là ma trận pháp Gauss-Jordan.Tuy nhiên, với các đơn vị. bài toán phi tuyến, hệ phương trình trên 3. n=n+1 không thể giải được. Khai triển F(X,α) 4. theo chuỗi Taylor, ta thu được biểu thức dưới dạng ma trận: 5. Tính (4) 6. Nếu n<2, đi tới bước 9 7. Nếu n<3, đi tới bước 8 Trong đó M là ma trận [m m] mà: 8. Nếu , trong đó thì tiếp tục (5 vòng lặp, nếu không, thoát ra ) khỏi vòng lặp. 9. Đặt ; Quay lại Và bước 2. 123
3. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 ắt đầu có n<2 không n<3 có Kết thúc Hình 1. Lưu đồ thuật toán 124
4. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 S d tr o Thủ tục làm khớp dữ liệu được thực hiện bởi hai chương trình con LSfit_NL (không có trọng số) và LSfit_NLW (có trọng số). Cú pháp khai báo như sau: LSfit_NL(matrix X, matrix Y, int par_num, matrix par) LSfit_NLW(matrix X, matrix Y, matrix W, int par_num, matrix par) Trong đó X, Y là hai ma trận tương ứng với bộ số liệu thực nghiệm (X,Y), W là ma trận trọng số, par_num là số tham số tự do của mô hình làm khớp, par là ma trận tương ứng với giá trị ban đầu của tham số. Mảng hai chiều hoặc một chiều có thể được chuyển thành ma trận (matrix) thông qua chương trình con array_to_matrix với cú pháp như sau: array_to_matrix((double *)array, int row, int col); array là mảng 1 chiều hoặc 2 chiều, row là số dòng, và col là số cột của ma trận tạo thành. Ví dụ, mảng hai chiều A[3][2] có thể được chuyển đổi thành ma trận MA[3][2] thông qua câu lệnh sau: MA = array_to_matrix((double*)A, 3, 2) Mô hình làm khớp được khai báo bên trong chương trình con uf double uf(double x, matrix par) { double result; result = par.E[0]*exp(x/par.E[1]); //par.E[i] là tham số tự do thứ i, x là biến. return result; } Đoạn chương trình thực hiện tác vụ làm khớp bộ số liệu 1 điểm theo mô hình f(x)=a*exp(x/b) với a, b là các tham số tự do được đưa ra dưới đây: double uf(double x, matrix par) #include #include "matrix.h" #include using namespace std; { double result; result = par.E[0]*exp(x/par.E[1]); // par.E[0]=a; par.E[1]=b return result; 125
5. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 } int main() { double A[15]={1,2,3,4,5,6,7,8,9,10,11,12,13,14,15}; double B[15]={20,39,66,113,180,300,497,816,1346,2230,3674,6050,9 976,16454,27122}; double C[15] = {0.222851, 0.160098, 0.122793, 0.094265, 0.074557 ,0.057718, 0.044859, 0.035006, 0.027256, 0.021178, 0.016497, 0.012857, 0.010012, 0.007796, 0.006072}; double para[2]={10,2}; matrix X = array_to_matrix(A,15,1); matrix Y = array_to_matrix(B,15,1); matrix W = array_to_matrix(C,15,1); matrix par = array_to_matrix(para,2,1); cout<<"No Weighted:"<<endl; Mprint(LSfit_NL(X,Y,2,par)); par=array_to_matrix(para,2,1); //khởi tạo lại tham số ban đầu cout<<"Weighted:"<<endl; Mprint(LSfit_NLW(X,Y,W,2,par)); return 0; } Mảng A, B, C lần lượt tương ứng với các ma trận X, Y, W. Để hiệu lực hóa chương trình, kết quả tính toán thực hiện bởi chương trình trên nhiều bộ số liệu khác nhau với các mô hình liệt kê dưới đây được so sánh với chương trình ROOT. Các mô hình làm khớp được thử nghiệm bao gồm: - Mô hình hàm lũy thừa cơ số tự nhiên: f(x)=a*exp(x/b); - Mô hình hàm gauss g(x)=A*exp(-(x- )2/2 ), với A, , là các tham số tự do. - Mô hình hàm gauss nằm trên một nền phông tương ứng với đa thức bậc 1: f(x) = g(x) + a1*x + a0 - Mô hình hai hàm gauss nằm chập lên nhau chồng trên một nền phông tương ứng với đa thức bậc 1: f(x) = g1(x) + g2(x)+ a1*x + a0 126
6. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 - Mô hình ba hàm gauss chập trên nền phông tương ứng với đa thức bậc 1: f(x)=f(x) = g1(x) + g2(x)+g3(x)+ a1*x + a0 3. H ệu tr t qu so s vớ ROOT Kết quả thu được bởi chương trình được so sánh với ROOT, phần mềm được sử dụng rộng rãi bởi nhiều phòng thí nghiệm trên thế giới. Kết quả so sánh với một số mô hình được trình bày trong Bảng 1. Bảng 1. So sánh giá trị tham số làm khớp của chương trình với ROOT Độ ệ C tr T m số Gí trị b đầu ROOT này tr (%) H m ũy t ừ số t ê : y= *exp(x/b) Không trọng số a 10 15,0015 15,0015 0 b 2 2,0000 2,0000 0 Có trọng số a 10 14,9828 15,0119 0,19 b 2 1,9997 2,0002 0,03 Hàm gauss: y = A*exp(-(x-μ)2/σ2) Không trọng số A 80 99,7945 99,7951 0,00 μ 52 50,0195 50,0195 0,00 σ 20 10,0041 10,0041 0,00 Có trọng số A 80 99,3681 99,9418 0,58 μ 52 50,1327 50,0060 0,25 σ 20 10,0439 9,9759 0,68 G uss + đ t ứ bậ 1: y = A*exp(-(x-μ)2/σ2)+ a1*x + a0 Không trọng số A 80 99,8306 99,8308 0,00 μ 52 50,0090 50,0090 0,00 σ 10 10,0198 10,0198 0,00 a1 2 2,0047 2,0047 0,00 a0 3 2,7612 2,7611 0,00 Có trọng số A 80 99,7179 99,9112 0,19 μ 52 50,0157 50,0015 0,03 σ 10 10,0227 10,0274 0,05 a1 2 2,0026 2,0076 0,25 a0 3 2,8558 2,5562 10,49 127
7. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 2 m uss ập + đ t ứ bậ 1: y = A*exp(-(x-μ)2/σ2)+ a1*x + a0 + A1*exp(-(x- μ1)2/σ12) Không trọng số A1 80 79,8317 79,8321 0,00 μ1 31 29,9800 29,9801 0,00 σ1 8 7,9845 7,9845 0,00 a1 2 0,9997 0,9997 0,00 a0 3 2,0129 2,0128 0,00 A2 90 99,9497 99,9500 μ2 51 49,9717 49,9717 σ2 10 12,0307 12,0307 Có trọng số A1 80 79,6992 79,8282 0,16 μ1 31 29,9754 29,9798 0,01 σ1 8 7,9896 7,9745 0,19 a1 2 0,9999 0,9997 0,03 a0 3 2,0053 2,0195 0,71 A2 90 99,8090 99,9412 0,13 μ2 51 49,9667 49,9676 0,00 σ2 10 12,0467 12,0371 0,08 3 m uss ập + đ t ứ bậ 1: A1*exp(-(x-μ1)2/σ12)+ a1*x + a0 + A2*exp(-(x- μ2)2/σ22)+ A3*exp(-(x-μ3)2/σ32) Không trọng số A1 75 80,1944 80,2094 0,02 μ1 25 31,0066 31,0096 0,01 σ1 10 7,9983 8,0012 0,04 a1 1 1,9997 1,9998 0,00 a0 2 2,9995 2,9983 0,04 A2 100 90,0007 89,9736 0,03 μ2 45 50,9908 50,9791 0,02 σ2 8 9,9805 9,9628 0,18 A3 210 200,2460 200,4460 0,10 μ3 65 60,9977 60,9968 0,00 σ3 8 4,9986 5,0010 0,05 Có trọng số A1 75 80,6545 80,3449 0,38 μ1 25 31,1255 31,0127 0,36 σ1 10 8,1015 8,0189 1,02 a1 1 2,0004 2,0020 0,08 a0 2 2,9873 2,8098 5,94 A2 100 88,8030 90,0044 1,35 128
8. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 μ2 45 50,5029 50,9623 0,91 σ2 8 9,2869 9,9319 6,95 A3 210 208,4220 200,7630 3,67 μ3 65 60,9426 60,9965 0,09 σ3 8 5,0936 5,0069 1,70 rất nhỏ. Ví dụ như trường hợp tham số 4. K t qu a thu được khi làm khớp với mô hình gauss trên nền đa thức bậc 1, độ lệch của Mô hình đầu tiên được sử dụng chương trình với ROOT là 10,4% để so sánh là mô hình hàm lũy thừa cơ (2,8585 so với 2,5562). Mặc dù độ lệch số tự nhiên. Mô hình hàm lũy thừa là cao, nhưng ảnh hưởng của tham số này dạng mô hình điển hình nhất khi tiến tới giá trị của hàm là rất nhỏ. hành thủ tục làm khớp phi tuyến, do tham số ảnh hưởng rất mạnh tới giá trị Kết quả có sự tương đồng cao của hàm. Chỉ một lượng nhỏ thay đổi giữa chương trình với ROOT khi áp trong tham số cũng khiến giá trị của dụng vào các mô hình chập gauss cho hàm thay đổi một lượng lớn.Kết quả thấy, chương trình hoàn toàn đáp ứng trong Bảng 1 cho thấy, khi làm khớp tốt bài toán tách đỉnh chập, vốn rất phổ không trọng số, chương trình hội tụ về biến khi phân tích phổ gamma. giá trị tham số hoàn toàn giống với 5. K t luận ROOT. Đối với quá trình làm khớp có trọng số, kết quả thulệch so với ROOT Chương trình làm khớp có kết một lượng nhỏ hơn .2 . quả có độ tương đồng cao với ROOT. Cấu trúc của chương trình đơn giản, Các mô hình gauss, gauss trên thuần túy chỉ sử dụng các thư viện có nền đa thức bậc một, chập 2 hàm gauss sẵn của C++, thuận tiện cho việc nhúng trên nền đa thức bậc 1, và chập 3 hàm vào các chương trình con khác. gauss trên nền đa thức bậc một đều cho kết quả tương đồng với ROOT. Độ Chương trình rất thích hợp để tích chênh lệch của giá trị tham số làm khớp hợp vào các chương trình phân tích phổ thu được bởi chương trình với giá trị thu tự thiết kế, qua đó giúp giảm chi phí mua được từ ROOT phần lớn đều nhỏ hơn phần mềm phân tích đắt tiền, với các tính 1%. Chỉ có một số ít trường hợp, giá trị năng ít hoặc không bao giờ được sử dụng. tham số làm khớp thu được bởi chương Ngoài ra, việc dễ dàng chỉnh sửa mã trình lệch so với ROOT cao hơn 1 . nguồn, giúp người dùng dễ dàng xây Tuy nhiên trong các trường hợp đó, các dựng các mô-đun chuyên biệt nhằm thực tham số có độ lệch cao là các tham số có hiện các tác vụ theo yêu cầu cụ thể một mức độ ảnh hưởng tới giá trị của hàm số cách thuận tiện và nhanh chóng. 129
9. TẠP CHÍ KHOA HỌC - ĐẠI HỌC ĐỒNG NAI, SỐ 03 - 2016 ISSN 2354-1482 TÀI LIỆU THAM KHẢO 1. [Online]. Available: 2. [Online]. Available: 3. [Online]. Available: 4. [Online]. Available: 5. [Online]. Available: 6. [Online]. Available: 7. Gill, P. R.; Murray, W.; and Wright, M. H. "The Levenberg-Marquardt Method." §4.7.3 in Practical Optimization. London: Academic Press, pp. 136-137, 1981 8. [Online]. Available: 9. Rao, C. R.; Toutenburg, H.; et al. (2008). Linear Models: Least Squares and Alternatives. Springer Series in Statistics (3rd ed.). Berlin: Springer. ISBN 978- 3-540-74226-5 DEVELOPMENT OF SUBROUTINE FOR DATA FITTING WITH VARIOUS MODELS ABSTRACT Currently, there are many computer programs, which allow users to fit experimental data to any mathematical models. However, these programs either do not give users their source codes (commercial software) or have complicated libraries (open source software). Consequently, using their libraries to form homemade software becomes a difficult task, and even impossible, in case of commercial software. This work presents a group of sub-programs, written in C++, which permit users to fit experimental data to any mathematical models, including weighted fit and non-weighted fit. The sub-programs are packaged in one file with only 438 code lines; hence, make it easy to develop programs based on these sub- programs. The quality of these sub-programs was proved by comparing with ROOT. Keywords:Fitting C++ code, Levenberg–Marquardt algorithm 130