Làm thế nào để bạn phù hợp với một gaussian trong python?

Ngày 19 tháng 11 năm 2018

Kiểm tra mã

Sự phong phú của phần mềm có sẵn để giúp bạn điều chỉnh các đỉnh vô tình làm phức tạp quá trình bằng cách chôn vùi các chức năng điều chỉnh toán học tương đối đơn giản dưới các lớp tính năng GUI. Như tôi hy vọng bạn đã thấy trong Phần I của loạt bài này, có thể khá mạnh mẽ để có thể phù hợp với tất cả dữ liệu của bạn từ đầu thay vì sử dụng phần mềm phù hợp nhất trong trường hợp này. Với bài đăng này, tôi muốn tiếp tục truyền cảm hứng cho bạn từ bỏ GUI và sử dụng python để xử lý dữ liệu của bạn bằng cách chỉ cho bạn cách khớp các đỉnh quang phổ với hình dạng đường và trích xuất nhiều thông tin để hỗ trợ phân tích của bạn

Tôi đã thêm sổ ghi chép mà tôi đã sử dụng để tạo bài đăng trên blog này, 181119_PeakFitting, vào kho lưu trữ GitHub của mình có thể tìm thấy tại đây. Cũng như bài đăng trước, bài đăng này được thiết kế để người đọc có thể theo dõi trong sổ ghi chép, và do đó, bài đăng này sẽ giải thích chức năng/nghĩa của từng ô thay vì cho bạn biết phải nhập gì cho từng ô

Đường nét Gaussian

Trước tiên, chúng tôi sẽ tập trung vào việc khớp các đường cong gaussian đơn và nhiều. Đầu tiên, tôi tạo một số dữ liệu gaussian giả để làm việc (xem sổ ghi chép và bài đăng trước)

Đường cong gaussian đơn.

Như bạn có thể thấy, điều này tạo ra một đỉnh duy nhất có dạng đường gaussian, với tâm, biên độ và chiều rộng cụ thể. Sau đó, chúng tôi muốn điều chỉnh đỉnh này thành một đường cong gaussian duy nhất để chúng tôi có thể trích xuất ba tham số này. Để làm như vậy, giống như với các đường cong tuyến tính hoặc hàm mũ, chúng tôi xác định hàm phù hợp mà chúng tôi sẽ đưa vào hàm scipy để khớp với dữ liệu giả.

def _1gaussian(x, amp1,cen1,sigma1):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen1)/sigma1)**2)))

Tôi đã xây dựng hàm phù hợp này bằng cách sử dụng phương trình cơ bản của phân phối gaussian. Sau đó, chúng tôi đưa hàm này vào một hàm scipy, cùng với dữ liệu trục x và trục y cũng như dự đoán của chúng tôi về các tham số phù hợp với hàm (tôi sử dụng các giá trị trung tâm, biên độ và sigma mà tôi đã sử dụng để tạo dữ liệu giả mạo

popt_gauss, pcov_gauss = scipy.optimize.curve_fit(_1gaussian, x_array, y_array_gauss, p0=[amp1, cen1, sigma1])
perr_gauss = np.sqrt(np.diag(pcov_gauss))

Sau đó, chúng tôi có thể in ra ba tham số phù hợp với các lỗi tương ứng của chúng

amplitude = 122.80 (+/-) 3.00
center = 49.90 (+/-) 0.33
sigma = 11.78 (+/-) 0.33

Và sau đó vẽ dữ liệu của chúng tôi cùng với sự phù hợp

Vừa đường cong gaussian đơn.

Điều chỉnh này thực hiện khá tốt công việc điều chỉnh dữ liệu gaussian giả. Bây giờ chúng ta có thể điều chỉnh thành công một đỉnh gaussian đơn được phân giải tốt, hãy làm việc với trường hợp phức tạp hơn khi chúng ta có một số đỉnh chồng chéo cần phải được kết hợp với nhau.

Điều cực kỳ hữu ích là có thể khớp một số đỉnh chồng chéo, bởi vì thông thường các đặc điểm quang phổ không được phân giải rõ ràng với nhau và thậm chí một phần đuôi của đường cong gaussian có thể làm lệch sự khớp của một đường cong khác nếu chúng chồng lên nhau. Chúng tôi sẽ xem xét một vùng dữ liệu chồng chéo mạnh mẽ

Hai đường cong gaussian chồng lên nhau.

Trong dữ liệu giả ở trên, bạn có thể thấy rằng có hai đỉnh gaussian có đặc điểm (i. e. tâm và chiều rộng của đỉnh) không thể phân biệt được với nhau. Trong nghiên cứu của mình, tôi thường gặp trường hợp này khi quan sát quang phổ Huỳnh quang của các loài phát ra ở bước sóng tương tự, và tôi quan tâm đến diện tích tương đối dưới hai đường cong; .

Giống như với sự phù hợp theo cấp số nhân mà chúng tôi đã nghiên cứu trước đây, để phù hợp với các đỉnh gaussian chồng chéo, chúng ta cần xác định một hàm cho tổng của hai gaussian

def _2gaussian(x, amp1,cen1,sigma1, amp2,cen2,sigma2):
    return amp1*(1/(sigma1*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen1)/sigma1)**2))) + \
            amp2*(1/(sigma2*(np.sqrt(2*np.pi))))*(np.exp((-1.0/2.0)*(((x_array-cen2)/sigma2)**2)))

Bạn có thể thấy rằng sự khác biệt duy nhất giữa _1gaussian và _2gaussian là cái sau là tổng của hai hàm gaussian và khớp với sáu tham số thay vì ba trong cái trước

Đưa cái này vào scipy hoạt động theo cùng một cách

popt_2gauss, pcov_2gauss = scipy.optimize.curve_fit(_2gaussian, x_array, y_array_2gauss, p0=[amp1, cen1, sigma1, amp2, cen2, sigma2])
perr_2gauss = np.sqrt(np.diag(pcov_2gauss))
pars_1 = popt_2gauss[0:3]
pars_2 = popt_2gauss[3:6]
gauss_peak_1 = _1gaussian(x_array, *pars_1)
gauss_peak_2 = _1gaussian(x_array, *pars_2)

Tuy nhiên, tôi đã thêm một vài dòng mã để xác định pars_1, pars_2 và gauss_peak_1, gauss_peak_2. Các biến phân tích cú pháp là các mảng chỉ chứa các tham số phù hợp cho các đỉnh thứ nhất và thứ hai tương ứng. Chúng tôi muốn xác định các biến này vì sau đó chúng tôi có thể đưa chúng vào hàm _1gaussian của mình để xây dựng dữ liệu cho các đỉnh độc lập

Nếu chúng ta vẽ biểu đồ dữ liệu hai gaussian giả của mình và khớp _2gaussian, chúng ta sẽ thấy rằng dữ liệu (các chấm đỏ) được theo dõi độc đáo bởi khớp (đường đứt nét màu đen)

Khớp hai gaussian.

Tuy nhiên, chúng tôi muốn có thể tự mình nhìn thấy các đỉnh sau khi chúng được tách ra khỏi nhau. Đây là nơi các biến gauss_peak_1 và _2 phát huy tác dụng.

Nếu chúng ta thêm các dòng mã sau vào ô biểu đồ của mình, chúng ta có thể tự vẽ hai đỉnh

ax1.plot(x_array, gauss_peak_1, "g")
ax1.fill_between(x_array, gauss_peak_1.min(), gauss_peak_1, facecolor="green", alpha=0.5)
ax1.plot(x_array, gauss_peak_2, "y")
ax1.fill_between(x_array, gauss_peak_2.min(), gauss_peak_2, facecolor="yellow", alpha=0.5)

Điều này sử dụng một chức năng gọi là “fill_between” từ thư viện matplotlib. Tôi thích sử dụng hàm này vì nó cho phép chúng ta thấy diện tích mà mỗi đỉnh chiếm dưới đường cong tổng

Các đường cong gaussian chồng chéo được giải mã.

Sau đó, chúng tôi có thể in các tham số phù hợp của hai đường cong, cũng như diện tích bên dưới các đường cong để trích xuất tất cả thông tin mong muốn nhằm phân tích hệ thống của chúng tôi.

-------------Peak 1-------------
amplitude = 120.09 (+/-) 3.54
center = 39.80 (+/-) 0.37
sigma = 11.39 (+/-) 0.41
area = 59.42
--------------------------------
-------------Peak 2-------------
amplitude = 78.51 (+/-) 2.73
center = 65.22 (+/-) 0.16
sigma = 5.15 (+/-) 0.16
area = 38.86
--------------------------------

Một cách phổ biến (quan trọng) để trực quan hóa dữ liệu của bạn là vẽ biểu đồ phần còn lại của dữ liệu cùng với sự phù hợp. Phần còn lại chỉ là những gì nó nghe như thế nào; . Nói cách khác, phần dữ liệu không phù hợp. Để tìm số dư, chúng ta có thể làm như sau

residual_2gauss = y_array_2gauss - (_2gaussian(x_array, *popt_2gauss))

Hơn nữa, chúng tôi lấy dữ liệu trục y (y_array_2gauss) và trừ đi độ khớp từ nó (_2gaussian(x_array, *popt_2gauss)) và gán giá trị này cho số dư_2gauss. Tôi muốn thêm phần này vào hình với các tín hiệu đã giải mã dưới dạng ô con bên dưới dữ liệu phù hợp ban đầu

Các đường cong gaussian chồng chéo được giải mã với phần dư.

Hình này là bức tranh hoàn chỉnh về kết quả của việc tách hai đường cong gaussian khỏi nhau. Tôi hy vọng rằng bạn có thể thấy cách bạn có thể dịch quy trình này thành khớp nhiều hơn hai đường cong cho quang phổ phức tạp hơn.

Bây giờ chúng ta đã có các đường hình dạng gaussian dưới vành đai của mình, tôi muốn nhanh chóng chỉ cho bạn cách khớp các đỉnh hình đường Lorentzian và Voigt. Để làm như vậy, tôi sẽ thực hiện theo cùng một quy trình chính xác, ngoại trừ tôi sẽ xác định các hàm khớp mới của loại Lorentzian và Voigt

Đường Lorentzian

Để cung cấp cho bạn thêm thực hành/ví dụ về khớp đỉnh, tôi sẽ minh họa cách khớp các đỉnh Lorentzian với ba đỉnh chồng lên nhau. Đây là một ví dụ về loại dữ liệu thu được từ quang phổ NMR, trong đó các đỉnh có hình dạng đường Lorentzian và thường có các bội số chồng lên nhau của các đỉnh

Các đường cong lorentz chồng lên nhau.

Hàm khớp mà chúng ta sử dụng ở đây là tính tổng của ba hàm Lorentzian.

________số 8

Khi chúng tôi khớp dữ liệu ba lorentzian với chức năng này, chúng tôi có thể in ra thông tin sau về ba đỉnh

-------------Peak 1-------------
amplitude = 50.26 (+/-) 1.07
center = 100.11 (+/-) 0.13
width = 5.84 (+/-) 0.18
area = 746.81
--------------------------------
-------------Peak 2-------------
amplitude = 100.45 (+/-) 0.80
center = 150.12 (+/-) 0.08
width = 10.60 (+/-) 0.13
area = 2659.01
--------------------------------
-------------Peak 3-------------
amplitude = 50.81 (+/-) 1.10
center = 200.11 (+/-) 0.12
width = 5.57 (+/-) 0.18
area = 720.14
--------------------------------

Sau đó, nếu chúng ta giải phần dư và vẽ đồ thị thông tin khớp tổng của mình, chúng ta có thể thấy rằng hàm khớp này thực hiện công việc khớp dữ liệu khá tốt

Giải chập các đường cong Lorentzian chồng lên nhau.

Như bạn có thể thấy, việc khớp các đỉnh hình dạng đường Lorentzian rất giống với các đỉnh gaussian, hãy lưu chức năng khớp. Cuối cùng, chúng ta sẽ xem xét cách điều chỉnh đỉnh hình dạng đường Voigt.

Không có gì đáng ngạc nhiên, việc điều chỉnh các đỉnh hình dạng đường Voigt tuân theo quy trình tương tự như điều chỉnh tất cả các loại dữ liệu khác mà chúng tôi đã xem xét cho đến nay. Trong trường hợp bạn chưa biết, các đường nét Voigt là sự kết hợp của các đường nét Lorentzian và Gaussian. Loại cực đại này phổ biến trong quang phổ nhiễu xạ tia X, trong đó chúng ta làm nhiễu xạ chùm tia X ra khỏi vật liệu rắn. Điều này có thể làm cho đỉnh mở rộng theo cả cơ chế Gaussian và Lorentzian

Dữ liệu hình dạng đường Voigt.

Bạn có thể thấy đỉnh nhọn hơn như thế nào, đây là một đặc điểm của đỉnh Lorentzian, trong khi khi bạn tiến gần đến đường cơ sở, đỉnh mở rộng ra, một đặc điểm của các đường cong Gaussian (i. e. phân phối chuẩn).

Do đó, hàm khớp mà chúng tôi sử dụng để khớp với đỉnh hình dạng đường Voigt là tổng trọng số của các hàm Gaussian và Lorentzian

popt_gauss, pcov_gauss = scipy.optimize.curve_fit(_1gaussian, x_array, y_array_gauss, p0=[amp1, cen1, sigma1])
perr_gauss = np.sqrt(np.diag(pcov_gauss))

0

Khi chúng tôi đưa điều này vào chức năng điều chỉnh scipy, như mong đợi, chúng tôi thấy rằng chức năng này phù hợp với dữ liệu khá tốt

Dữ liệu hình dạng đường Voigt phù hợp với hàm Voigt.

Một điều cuối cùng cần chỉ ra về chức năng điều chỉnh Voigt này là bởi vì chúng ta có các tham số trọng số Gaussian và Lorentzian (i. e. biên độ), chúng ta có thể sử dụng chúng để phân tích mức độ mở rộng cực đại từ các cơ chế mở rộng tương ứng. Điều đó khá tuyệt.

popt_gauss, pcov_gauss = scipy.optimize.curve_fit(_1gaussian, x_array, y_array_gauss, p0=[amp1, cen1, sigma1])
perr_gauss = np.sqrt(np.diag(pcov_gauss))

1

Nếu bạn đã tiến xa đến mức này trong loạt bài phù hợp, xin chúc mừng. Giờ đây, bạn đã có tất cả các kỹ năng cơ bản (và cả một số sổ tay tuyệt vời nữa) để phù hợp với dữ liệu thử nghiệm của mình. Tôi sẽ tiếp tục loạt bài này trong tương lai để đi vào chi tiết cụ thể về việc khớp các loại dữ liệu quang phổ khác nhau, vì vậy hãy chú ý theo dõi nếu bạn đang sử dụng quang phổ Huỳnh quang, NMR hoặc XRD, vì tôi sẽ đề cập đến việc khớp loại dữ liệu này trong các bài đăng trong tương lai

---

---

Linkedin Github Instagram Youtube

Tất cả những suy nghĩ và ý kiến ​​​​là của riêng tôi và không phản ánh những suy nghĩ và ý kiến ​​​​của tổ chức của tôi

Phù hợp Gaussian là gì?

Curve_fit hoạt động như thế nào với Python?