戰(zhàn)：從原理到解決重疊峰分解難題)
1. 項(xiàng)目概述多峰高斯擬合的挑戰(zhàn)與價(jià)值在信號(hào)處理、光譜分析、色譜分離、生物醫(yī)學(xué)成像乃至金融數(shù)據(jù)分析中我們常常會(huì)遇到一種典型的數(shù)據(jù)形態(tài)一個(gè)看似復(fù)雜的波形實(shí)際上是由多個(gè)獨(dú)立的、相互重疊的“峰”疊加而成。比如一張質(zhì)譜圖上的多個(gè)離子峰一條光譜中不同元素的特征譜線(xiàn)或者心電圖里相鄰的P波、QRS波和T波。直接觀察這些混合在一起的峰我們很難精確地知道每個(gè)峰的中心位置、高度和寬度而這些參數(shù)恰恰是定量分析的核心。這時(shí)多峰高斯擬合就成了從混沌中提取秩序的“數(shù)學(xué)手術(shù)刀”。高斯函數(shù)或者說(shuō)正態(tài)分布曲線(xiàn)因其完美的鐘形對(duì)稱(chēng)性和良好的數(shù)學(xué)性質(zhì)成為描述這些獨(dú)立峰最常用的模型。單峰擬合相對(duì)簡(jiǎn)單但當(dāng)多個(gè)高斯峰擠在一起尤其是高度不同、寬度不一、甚至基線(xiàn)還有傾斜或偏移時(shí)問(wèn)題就變得棘手了。手動(dòng)“猜”參數(shù)幾乎不可能而簡(jiǎn)單的自動(dòng)擬合算法很容易陷入局部最優(yōu)解給出完全不合理的結(jié)果比如把兩個(gè)峰擬合成了一個(gè)寬峰或者擬合出的峰跑到數(shù)據(jù)范圍之外去了。我這次要分享的就是在MATLAB環(huán)境下成功實(shí)現(xiàn)三個(gè)重疊高斯峰的精確擬合全過(guò)程。這不僅僅是調(diào)用一個(gè)fit函數(shù)那么簡(jiǎn)單它涉及對(duì)數(shù)據(jù)本質(zhì)的理解、初始參數(shù)的巧妙估計(jì)、擬合算法的選擇以及大量“踩坑”后總結(jié)出的調(diào)試技巧。無(wú)論你是分析實(shí)驗(yàn)數(shù)據(jù)的研究生還是處理監(jiān)測(cè)數(shù)據(jù)的工程師這套從數(shù)據(jù)預(yù)處理、模型構(gòu)建、參數(shù)初始化到結(jié)果驗(yàn)證的完整流程都能讓你在面對(duì)復(fù)雜多峰數(shù)據(jù)時(shí)心里更有底。2. 核心思路與模型構(gòu)建理解“擬合”在做什么在動(dòng)手寫(xiě)代碼之前我們必須徹底想清楚多峰高斯擬合我們到底在求什么這決定了我們整個(gè)方案的架構(gòu)。2.1 數(shù)學(xué)模型三個(gè)高斯峰的疊加我們的目標(biāo)模型是三個(gè)高斯函數(shù)的線(xiàn)性疊加再加上一個(gè)可能存在的基線(xiàn)Baseline。一個(gè)標(biāo)準(zhǔn)的高斯函數(shù)公式如下y A * exp(-(x - μ)^2 / (2 * σ^2))其中A 峰高Amplitude。決定了峰的最大值。μ 峰位Mean/Center。決定了峰在x軸上的中心位置。σ 標(biāo)準(zhǔn)差Standard Deviation。決定了峰的寬度。半高全寬FWHM與σ的關(guān)系為FWHM 2√(2ln2) * σ ≈ 2.355 * σ。對(duì)于三個(gè)峰我們的總模型就是y_total Baseline Peak1 Peak2 Peak3即y_total (b0 b1*x) A1*exp(-(x-μ1)^2/(2*σ1^2)) A2*exp(-(x-μ2)^2/(2*σ2^2)) A3*exp(-(x-μ3)^2/(2*σ3^2))這里我特意將基線(xiàn)設(shè)為一次線(xiàn)性項(xiàng)b0 b1*x而不是一個(gè)常數(shù)。這是因?yàn)樵趯?shí)際數(shù)據(jù)中特別是光譜或色譜數(shù)據(jù)由于儀器背景或漂移基線(xiàn)傾斜非常常見(jiàn)。忽略它會(huì)導(dǎo)致峰高和峰面積的估計(jì)產(chǎn)生系統(tǒng)誤差。注意是否包含基線(xiàn)、基線(xiàn)是常數(shù)b0還是一次項(xiàng)b0b1*x甚至二次項(xiàng)需要根據(jù)你的數(shù)據(jù)實(shí)際情況判斷。一個(gè)簡(jiǎn)單的判斷方法是觀察數(shù)據(jù)中“無(wú)峰”區(qū)域的趨勢(shì)。如果拿不準(zhǔn)從簡(jiǎn)單模型常數(shù)基線(xiàn)開(kāi)始嘗試如果擬合殘差數(shù)據(jù)點(diǎn)與擬合曲線(xiàn)的差值呈現(xiàn)明顯的趨勢(shì)性分布則說(shuō)明需要更復(fù)雜的基線(xiàn)模型。2.2 擬合的本質(zhì)非線(xiàn)性最小二乘優(yōu)化擬合就是尋找一組模型參數(shù)上面提到的A1, μ1, σ1, A2, μ2, σ2, A3, μ3, σ3, b0, b1使得模型計(jì)算出的曲線(xiàn)y_total與實(shí)測(cè)數(shù)據(jù)點(diǎn)y_data之間的總體差異最小。這個(gè)差異通常用殘差平方和RSS來(lái)衡量RSS Σ(y_data_i - y_total_i)^2。擬合過(guò)程就是一個(gè)不斷調(diào)整這11個(gè)參數(shù)讓RSS達(dá)到最小的優(yōu)化過(guò)程。由于高斯函數(shù)是非線(xiàn)性的這是一個(gè)“非線(xiàn)性最小二乘”問(wèn)題。MATLAB的lsqcurvefit或曲線(xiàn)擬合工具箱的fit函數(shù)內(nèi)部使用的就是諸如“Levenberg-Marquardt”之類(lèi)的算法來(lái)解決這個(gè)問(wèn)題。這類(lèi)算法非常強(qiáng)大但有一個(gè)致命弱點(diǎn)高度依賴(lài)初始參數(shù)猜測(cè)。如果初始值離真實(shí)值太遠(yuǎn)算法極易收斂到錯(cuò)誤的局部最優(yōu)解而不是全局最優(yōu)解。因此整個(gè)多峰擬合成功的關(guān)鍵一半在于構(gòu)建正確的模型另一半就在于如何為這11個(gè)參數(shù)提供一個(gè)“聰明”的初始估計(jì)。接下來(lái)我們就進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)。3. 數(shù)據(jù)準(zhǔn)備與初始參數(shù)估計(jì)為成功擬合奠基假設(shè)我們有一組實(shí)測(cè)數(shù)據(jù)x是自變量如波長(zhǎng)、時(shí)間、質(zhì)量數(shù)y是因變量如強(qiáng)度、吸光度、響應(yīng)值。數(shù)據(jù)已經(jīng)以數(shù)組形式存在于MATLAB工作區(qū)。3.1 數(shù)據(jù)可視化與初步觀察第一步永遠(yuǎn)是把數(shù)據(jù)畫(huà)出來(lái)用肉眼進(jìn)行第一次“診斷”。figure; plot(x, y, b.-, LineWidth, 1, MarkerSize, 10); xlabel(X (e.g., Wavelength)); ylabel(Y (e.g., Intensity)); title(Raw Data - Visual Inspection for Peaks); grid on;仔細(xì)觀察圖形識(shí)別峰的數(shù)量目標(biāo)是3個(gè)但你要確認(rèn)數(shù)據(jù)中是否明顯有3個(gè)凸起。有時(shí)噪聲或畸變會(huì)產(chǎn)生“假峰”。估計(jì)峰位μ用鼠標(biāo)光標(biāo)大致讀取三個(gè)峰頂對(duì)應(yīng)的x坐標(biāo)。記下它們例如mu1_guess, mu2_guess, mu3_guess。估計(jì)峰高A大致估計(jì)每個(gè)峰頂?shù)膟值。注意由于峰重疊這個(gè)值會(huì)低于該峰獨(dú)立存在時(shí)的高度?？梢韵扔梅逯禍p去附近“谷底”的y值來(lái)粗略估計(jì)。觀察基線(xiàn)看看數(shù)據(jù)最左邊和最右邊的點(diǎn)以及峰谷之間的區(qū)域連線(xiàn)趨勢(shì)是水平的還是傾斜的這決定了基線(xiàn)模型。3.2. 自動(dòng)化初始估計(jì)技巧對(duì)于更復(fù)雜或批量的數(shù)據(jù)我們可以借助MATLAB內(nèi)置函數(shù)進(jìn)行輔助估計(jì)這比肉眼估計(jì)更穩(wěn)健。1. 尋找峰值點(diǎn)估計(jì) μ 和 A使用findpeaks函數(shù)。這個(gè)函數(shù)能幫你找到局部極大值點(diǎn)并忽略一些小噪聲峰。[pks, locs, widths, proms] findpeaks(y, x, ... MinPeakProminence, max(y)*0.05, ... % 設(shè)置最小峰突出度過(guò)濾噪聲 MinPeakDistance, range(x)*0.05); % 設(shè)置最小峰間距避免識(shí)別到同一個(gè)峰pks是找到的峰值高度A的初始估計(jì)。locs是對(duì)應(yīng)的x位置μ的初始估計(jì)。widths和proms可以輔助估計(jì)σ。 檢查找到的峰數(shù)量是否大于等于3。如果太多可以調(diào)整MinPeakProminence和MinPeakDistance參數(shù)如果太少則調(diào)低這些閾值。2. 估計(jì)峰寬σ高斯峰的寬度σ可以通過(guò)多種方式估計(jì)半高寬法對(duì)于一個(gè)孤立的峰找到峰高一半處的兩個(gè)點(diǎn)其x坐標(biāo)差值即為FWHM然后除以2.355得到σ。對(duì)于重疊峰這很難自動(dòng)完成。利用findpeaks的輸出findpeaks函數(shù)返回的widths是每個(gè)峰在半高處的寬度這近似就是FWHM。我們可以用sigma_guess widths / 2.355。經(jīng)驗(yàn)公式如果數(shù)據(jù)點(diǎn)間隔均勻可以觀察峰上升沿/下降沿的陡峭程度。一個(gè)粗略的起點(diǎn)是設(shè)σ的初始值為相鄰峰間距的1/5到1/10。3. 估計(jì)基線(xiàn)參數(shù)b0, b1一個(gè)簡(jiǎn)單有效的方法是對(duì)數(shù)據(jù)兩端例如前5%和后5%的數(shù)據(jù)點(diǎn)進(jìn)行線(xiàn)性擬合得到的截距和斜率作為b0和b1的初始值。n length(x); indices [1:round(n*0.05), round(n*0.95):n]; % 取頭尾5%的索引 p polyfit(x(indices), y(indices), 1); % 一階多項(xiàng)式擬合 b0_guess p(2); b1_guess p(1);假設(shè)通過(guò)以上方法我們得到了如下初始猜測(cè)% 峰1 (最左側(cè)) A1_guess pks(1); mu1_guess locs(1); sigma1_guess widths(1)/2.355; % 峰2 (中間) A2_guess pks(2); mu2_guess locs(2); sigma2_guess widths(2)/2.355; % 峰3 (最右側(cè)) A3_guess pks(3); mu3_guess locs(3); sigma3_guess widths(3)/2.355; % 基線(xiàn) % b0_guess, b1_guess 來(lái)自polyfit將這些初始值組合成一個(gè)向量initial_guess [A1_guess, mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess, A3_guess, mu3_guess, sigma3_guess, b0_guess, b1_guess]。4. 擬合實(shí)現(xiàn)兩種主流方法與詳細(xì)步驟有了模型和初始參數(shù)我們就可以開(kāi)始擬合了。這里介紹兩種最常用的方法使用lsqcurvefit優(yōu)化函數(shù)和使用曲線(xiàn)擬合工具箱的fit函數(shù)。lsqcurvefit更底層、靈活適合集成到腳本中fit函數(shù)更直觀、快捷適合交互式分析。4.1 方法一使用lsqcurvefit進(jìn)行擬合lsqcurvefit是優(yōu)化工具箱中的函數(shù)它直接處理非線(xiàn)性最小二乘問(wèn)題。第一步定義模型函數(shù)在MATLAB中創(chuàng)建一個(gè)函數(shù)文件例如multiGauss.m或者使用匿名函數(shù)。這里用匿名函數(shù)示例% 定義三峰高斯帶線(xiàn)性基線(xiàn)的模型函數(shù) % params: [A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1] multiGaussModel (params, x) ... params(1) * exp(-(x - params(2)).^2 / (2 * params(3)^2)) ... % 峰1 params(4) * exp(-(x - params(41)).^2 / (2 * params(42)^2)) ... % 峰2 (注意索引) params(7) * exp(-(x - params(71)).^2 / (2 * params(72)^2)) ... % 峰3 params(10) params(11) * x; % 線(xiàn)性基線(xiàn) b0 b1*x注意索引的對(duì)應(yīng)關(guān)系。為了清晰也可以將參數(shù)解包multiGaussModel (p, x) ... p(1)*exp(-(x-p(2)).^2/(2*p(3)^2)) ... p(4)*exp(-(x-p(5)).^2/(2*p(6)^2)) ... p(7)*exp(-(x-p(8)).^2/(2*p(9)^2)) ... p(10) p(11)*x;第二步設(shè)置邊界約束關(guān)鍵步驟這是避免擬合出荒謬結(jié)果如負(fù)的峰寬、峰位跑到天涯海角的關(guān)鍵。我們需要為每個(gè)參數(shù)設(shè)置合理的上下界lb和ub。% 基于初始猜測(cè)設(shè)置邊界 % 順序: [A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1] % 下界 (Lower Bounds) lb [0, min(x), 0, ... % 峰1: 振幅0, 峰位在數(shù)據(jù)范圍內(nèi) 寬度0 0, min(x), 0, ... % 峰2 0, min(x), 0, ... % 峰3 -inf, -inf]; % 基線(xiàn)參數(shù)可以為任意值 % 上界 (Upper Bounds) ub [inf, max(x), range(x)/2, ... % 峰1: 振幅無(wú)上限峰位在數(shù)據(jù)范圍內(nèi)寬度小于數(shù)據(jù)范圍一半合理假設(shè) inf, max(x), range(x)/2, ... % 峰2 inf, max(x), range(x)/2, ... % 峰3 inf, inf]; % 基線(xiàn)參數(shù)無(wú)限制 % 可以更精細(xì)地約束例如讓峰位按順序排列避免擬合時(shí)峰位互換 % lb(2) lb(5) lb(8) 且 ub(2) ub(5) ub(8) 的邏輯更復(fù)雜通?？亢玫某跏贾当苊?。第三步執(zhí)行擬合% 設(shè)置優(yōu)化選項(xiàng)提高顯示細(xì)節(jié) options optimoptions(lsqcurvefit, Display, iter, Algorithm, trust-region-reflective); % ‘levenberg-marquardt’算法不支持邊界這里用‘trust-region-reflective’ % 執(zhí)行擬合 [params_fitted, resnorm, residual, exitflag, output] ... lsqcurvefit(multiGaussModel, initial_guess, x, y, lb, ub, options); disp(擬合參數(shù) (A, mu, sigma, b0, b1):); disp(params_fitted);params_fitted就是擬合得到的最優(yōu)參數(shù)向量。4.2 方法二使用曲線(xiàn)擬合工具箱fit函數(shù)fit函數(shù)語(yǔ)法更貼近“擬合”這個(gè)概念并且能自動(dòng)生成豐富的統(tǒng)計(jì)信息和繪圖。第一步定義擬合類(lèi)型和選項(xiàng)% 使用 fittype 定義模型coefficients 指定參數(shù)名稱(chēng) ft fittype(A1*exp(-(x-mu1)^2/(2*sigma1^2)) A2*exp(-(x-mu2)^2/(2*sigma2^2)) A3*exp(-(x-mu3)^2/(2*sigma3^2)) b0 b1*x, ... independent, x, ... dependent, y, ... coefficients, {A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1}); % 設(shè)置擬合選項(xiàng)包括初始值和邊界 opts fitoptions(ft); opts.StartPoint initial_guess; % 傳入我們之前準(zhǔn)備好的初始猜測(cè)向量 opts.Lower lb; % 下界 opts.Upper ub; % 上界 opts.Display Iter; % 顯示迭代過(guò)程 % opts.Robust LAR; % 如果數(shù)據(jù)有異常點(diǎn)可以嘗試穩(wěn)健擬合第二步執(zhí)行擬合并繪圖% 執(zhí)行擬合 [fitresult, gof] fit(x, y, ft, opts); % 顯示擬合結(jié)果 disp(fitresult); disp(gof); % 輸出擬合優(yōu)度統(tǒng)計(jì)量如 R-square, RMSE % 繪制擬合結(jié)果 figure; plot(fitresult, x, y); legend(原始數(shù)據(jù), 擬合曲線(xiàn), Location, Best); xlabel(X); ylabel(Y); title(三峰高斯擬合結(jié)果);fitresult是一個(gè)包含所有擬合參數(shù)的對(duì)象可以通過(guò)fitresult.A1,fitresult.mu1等方式訪問(wèn)。gof包含了擬合優(yōu)度的信息如決定系數(shù)rsquare越接近1說(shuō)明擬合越好。5. 結(jié)果評(píng)估、可視化與問(wèn)題排查擬合完成并不意味著結(jié)束我們必須嚴(yán)格評(píng)估擬合質(zhì)量并診斷可能的問(wèn)題。5.1 可視化評(píng)估四象限診斷圖一張好的診斷圖勝過(guò)千言萬(wàn)語(yǔ)。我習(xí)慣同時(shí)繪制四個(gè)子圖figure(Position, [100, 100, 1200, 800]); % 子圖1原始數(shù)據(jù) vs. 擬合曲線(xiàn) subplot(2,2,1); plot(x, y, b., MarkerSize, 8); hold on; x_fine linspace(min(x), max(x), 1000); % 生成更密的點(diǎn)用于繪制光滑曲線(xiàn) y_fitted multiGaussModel(params_fitted, x_fine); % 或用 fitresult(x_fine) plot(x_fine, y_fitted, r-, LineWidth, 2); legend(Data, Fitted Curve, Location, Best); title(Fit Overview); grid on; % 子圖2殘差圖 (Residuals) subplot(2,2,2); residuals y - multiGaussModel(params_fitted, x); % 計(jì)算殘差 plot(x, residuals, k^, MarkerSize, 5, MarkerFaceColor, k); hold on; plot([min(x), max(x)], [0,0], r--); % 零參考線(xiàn) xlabel(X); ylabel(Residual); title(Residual Plot); grid on; % 好的擬合殘差應(yīng)隨機(jī)分布在零線(xiàn)附近無(wú)趨勢(shì)性。 % 子圖3殘差直方圖 subplot(2,2,3); histogram(residuals, 20, Normalization, probability, FaceColor, [0.5, 0.5, 0.5]); xlabel(Residual); ylabel(Probability); title(Residual Distribution); grid on; % 理想情況應(yīng)接近均值為0的正態(tài)分布。 % 子圖4分峰顯示 (Peak Decomposition) subplot(2,2,4); plot(x_fine, y_fitted, k-, LineWidth, 1.5); hold on; % 計(jì)算并繪制每個(gè)單獨(dú)的峰 peak1 params_fitted(1) * exp(-(x_fine - params_fitted(2)).^2 / (2 * params_fitted(3)^2)); peak2 params_fitted(4) * exp(-(x_fine - params_fitted(5)).^2 / (2 * params_fitted(6)^2)); peak3 params_fitted(7) * exp(-(x_fine - params_fitted(8)).^2 / (2 * params_fitted(9)^2)); baseline params_fitted(10) params_fitted(11) * x_fine; plot(x_fine, peak1, g--, LineWidth, 1); plot(x_fine, peak2, b--, LineWidth, 1); plot(x_fine, peak3, m--, LineWidth, 1); plot(x_fine, baseline, c:, LineWidth, 1); legend(Total Fit, Peak 1, Peak 2, Peak 3, Baseline, Location, Best); title(Peak Decomposition); grid on;通過(guò)這四張圖你可以一目了然地判斷總覽圖擬合曲線(xiàn)是否完美貼合數(shù)據(jù)點(diǎn)殘差圖殘差是否隨機(jī)、無(wú)規(guī)律如果呈現(xiàn)“U”型或“∩”型說(shuō)明模型選擇不當(dāng)如基線(xiàn)模型不對(duì)。殘差分布是否近似正態(tài)嚴(yán)重偏離可能暗示有異常點(diǎn)或模型系統(tǒng)誤差。分峰圖分解出的單個(gè)峰是否合理有沒(méi)有出現(xiàn)負(fù)峰圖形上表現(xiàn)為向下凸峰位是否與預(yù)期相符5.2 定量評(píng)估指標(biāo)除了看圖還要看數(shù)決定系數(shù) R2gof.rsquare。大于0.99通常說(shuō)明擬合很好但要注意對(duì)于非常尖銳的峰即使R2很高峰面積也可能不準(zhǔn)。殘差平方和 (RSS) 或 均方根誤差 (RMSE)sqrt(gof.sse / length(x))。越小越好但要在不同數(shù)據(jù)集間比較才有意義。參數(shù)置信區(qū)間使用confint(fitresult)可以計(jì)算參數(shù)的95%置信區(qū)間。如果某個(gè)參數(shù)的置信區(qū)間非常寬例如包含0或負(fù)值說(shuō)明該參數(shù)不可靠可能是數(shù)據(jù)信息不足或者該參數(shù)與其他參數(shù)強(qiáng)相關(guān)共線(xiàn)性。5.3 常見(jiàn)問(wèn)題與排查技巧實(shí)錄在實(shí)際操作中你幾乎一定會(huì)遇到下面這些問(wèn)題。以下是我的排查清單問(wèn)題1擬合失敗提示“未收斂”或“達(dá)到最大迭代次數(shù)”。原因初始值太差或者邊界設(shè)置不合理導(dǎo)致優(yōu)化算法找不到下降方向。解決放松邊界先將所有邊界設(shè)得非常寬如lb -inf(1,11); ub inf(1,11);只保留sigma0這樣的物理約束。如果能擬合再逐步收緊邊界。改進(jìn)初始值回到第3步用更穩(wěn)健的方法如對(duì)數(shù)據(jù)平滑后再找峰估計(jì)初始值?？梢試L試手動(dòng)在圖上選點(diǎn)。分步擬合先擬合一個(gè)峰固定其參數(shù)再加入第二個(gè)峰擬合以此類(lèi)推。這能有效降低優(yōu)化難度。換用算法lsqcurvefit可以嘗試‘levenberg-marquardt’算法但不支持邊界。fit函數(shù)可以嘗試‘Robust’選項(xiàng)。問(wèn)題2擬合結(jié)果中某個(gè)峰的振幅是負(fù)值或者峰寬極大/極小。原因典型的局部最優(yōu)解或者模型過(guò)于復(fù)雜過(guò)擬合。解決施加物理約束強(qiáng)制振幅A大于0峰寬σ在一個(gè)合理范圍內(nèi)如數(shù)據(jù)范圍的1/100到1/2。簡(jiǎn)化模型檢查是否真的需要三個(gè)峰也許兩個(gè)峰加一個(gè)更復(fù)雜的基線(xiàn)模型就夠了?；蛘邍L試固定其中一兩個(gè)你認(rèn)為最確定的參數(shù)如已知某個(gè)峰位。檢查數(shù)據(jù)質(zhì)量數(shù)據(jù)噪聲是否太大考慮先對(duì)原始數(shù)據(jù)進(jìn)行平滑處理如Savitzky-Golay濾波但注意平滑可能扭曲峰形。問(wèn)題3殘差圖顯示出明顯的規(guī)律性如彎曲趨勢(shì)。原因模型不足以描述數(shù)據(jù)。通常是基線(xiàn)模型不合適。解決升級(jí)基線(xiàn)模型將常數(shù)基線(xiàn)b0改為線(xiàn)性b0b1*x甚至二次b0b1*xb2*x^2。檢查峰函數(shù)模型數(shù)據(jù)峰形是否不對(duì)稱(chēng)高斯函數(shù)是對(duì)稱(chēng)的。如果峰有明顯拖尾可能需要考慮洛倫茲Lorentzian函數(shù)或二者的混合Voigt profile。此時(shí)模型應(yīng)改為A / (1 ((x-μ)/σ)^2)。問(wèn)題4兩個(gè)峰的參數(shù)擬合后幾乎一樣或者峰位互換了。原因初始值中兩個(gè)峰的估計(jì)位置太接近或者算法在迭代中發(fā)生了“跳變”。解決嚴(yán)格約束峰位順序在lsqcurvefit中可以通過(guò)設(shè)置非線(xiàn)性約束來(lái)實(shí)現(xiàn)但這比較復(fù)雜。更實(shí)用的方法是在初始值中明確指定mu1_guess mu2_guess mu3_guess并設(shè)置不重疊的邊界如[mu1_guess-Δ, mu1_guessΔ]。使用“鎖定”策略先擬合最左側(cè)和最右側(cè)的兩個(gè)峰將它們的參數(shù)固定再擬合中間的峰。問(wèn)題5擬合速度很慢尤其是數(shù)據(jù)點(diǎn)很多的時(shí)候。原因每次迭代都要計(jì)算整個(gè)模型函數(shù)數(shù)據(jù)點(diǎn)多則計(jì)算量大。解決數(shù)據(jù)降采樣在保持峰形特征的前提下對(duì)數(shù)據(jù)進(jìn)行適當(dāng)降采樣。提供解析雅可比矩陣對(duì)于lsqcurvefit可以提供一個(gè)函數(shù)來(lái)計(jì)算模型關(guān)于各個(gè)參數(shù)的導(dǎo)數(shù)雅可比矩陣這能極大加速收斂。但對(duì)于高斯模型手動(dòng)推導(dǎo)并編寫(xiě)雅可比矩陣比較繁瑣除非對(duì)性能有極致要求否則通常不需要。6. 進(jìn)階技巧與擴(kuò)展應(yīng)用當(dāng)你掌握了三峰擬合后可以嘗試以下進(jìn)階操作讓分析更上一層樓。6.1 自動(dòng)化與批處理如果你有成百上千條光譜需要分析手動(dòng)操作是不可行的。你需要將上述流程封裝成函數(shù)。function [fittedParams, gofStats, fitResult] fitThreeGaussPeaks(xData, yData, initialGuess, lowerBounds, upperBounds) % 封裝三峰高斯擬合流程 % 輸入xData, yData, 初始猜測(cè)下界上界 % 輸出擬合參數(shù)統(tǒng)計(jì)量fit結(jié)果對(duì)象 ft fittype(...); % 定義模型 opts fitoptions(...); % 設(shè)置選項(xiàng) % ... [填充具體代碼] [fitResult, gofStats] fit(xData, yData, ft, opts); fittedParams coeffvalues(fitResult); end然后在一個(gè)循環(huán)中調(diào)用這個(gè)函數(shù)處理每個(gè)數(shù)據(jù)文件。關(guān)鍵難點(diǎn)在于自動(dòng)生成可靠的初始猜測(cè)。你可以開(kāi)發(fā)一個(gè)穩(wěn)健的峰值檢測(cè)和基線(xiàn)估計(jì)子函數(shù)作為fitThreeGaussPeaks的前置步驟。6.2 從擬合參數(shù)到物理量峰面積計(jì)算在很多應(yīng)用中峰高A受儀器條件影響大而峰面積Area更能代表物質(zhì)的量。對(duì)于高斯峰面積S A * σ * √(2π)。A1 params_fitted(1); sigma1 params_fitted(3); area1 A1 * sigma1 * sqrt(2*pi); % 同理計(jì)算 area2, area3如果存在線(xiàn)性基線(xiàn)上述公式計(jì)算的是峰相對(duì)于基線(xiàn)的凈面積這是正確的。6.3 模型選擇高斯 vs. 洛倫茲 vs. 沃伊特不是所有的峰都是完美的高斯形。高斯峰源于多普勒增寬、某些色譜過(guò)程。峰形較“瘦”衰減更快。洛倫茲峰源于自然增寬、某些共振現(xiàn)象。峰形較“胖”有更長(zhǎng)的拖尾。沃伊特峰高斯和洛倫茲的卷積能描述更復(fù)雜的增寬機(jī)制。在MATLAB中只需修改fittype中的表達(dá)式即可切換模型洛倫茲‘A / (1 ((x-mu)/sigma)^2)’沃伊特需要自定義函數(shù)或使用Faddeeva函數(shù)近似較為復(fù)雜。選擇模型的依據(jù)一是物理過(guò)程的先驗(yàn)知識(shí)二是看哪種模型的殘差更小、更隨機(jī)。可以都試一下用gof.rsquare和殘差圖來(lái)輔助判斷。6.4 不確定性分析與誤差傳遞擬合出的參數(shù)是有不確定性的置信區(qū)間。當(dāng)我們用這些參數(shù)計(jì)算衍生量如峰面積、峰位差時(shí)誤差也會(huì)傳遞。 MATLAB的曲線(xiàn)擬合工具箱可以計(jì)算參數(shù)的雅可比矩陣和協(xié)方差矩陣。對(duì)于簡(jiǎn)單的線(xiàn)性函數(shù)如面積計(jì)算誤差傳遞可以用公式近似。但對(duì)于復(fù)雜情況推薦使用蒙特卡洛模擬假設(shè)擬合參數(shù)服從以最佳估計(jì)值為均值、以標(biāo)準(zhǔn)誤差為方差的多維正態(tài)分布。從這個(gè)分布中隨機(jī)抽取大量如10000組參數(shù)集。對(duì)每組參數(shù)計(jì)算你關(guān)心的衍生量如面積。衍生量結(jié)果的分布如直方圖就給出了該量的估計(jì)值及其不確定性如95%置信區(qū)間。這個(gè)過(guò)程在MATLAB中實(shí)現(xiàn)起來(lái)需要一些編程但它提供了最全面的不確定性評(píng)估。成功擬合三個(gè)重疊高斯峰是一個(gè)從理論到實(shí)踐再到經(jīng)驗(yàn)積累的完整過(guò)程。它考驗(yàn)的不僅僅是對(duì)MATLAB函數(shù)的熟悉程度更是對(duì)數(shù)據(jù)、模型和優(yōu)化算法的綜合理解。最深刻的體會(huì)是沒(méi)有“一鍵萬(wàn)能”的擬合。初始值的精心設(shè)置、物理約束的合理施加、以及基于殘差圖的模型診斷這些手動(dòng)步驟的重要性往往超過(guò)選擇哪個(gè)具體的擬合函數(shù)。每次擬合都是一次與數(shù)據(jù)的對(duì)話(huà)你需要不斷提出問(wèn)題模型對(duì)嗎初始值好嗎并根據(jù)數(shù)據(jù)的“回答”殘差圖、參數(shù)值來(lái)調(diào)整策略。當(dāng)你看到分解出的三個(gè)光滑鐘形曲線(xiàn)完美地拼合成原始數(shù)據(jù)時(shí)那種從雜亂中提煉出清晰信息的成就感正是數(shù)據(jù)分析工作最大的樂(lè)趣所在。 MATLAB多峰高斯擬合實(shí)戰(zhàn)：從原理到解決重疊峰分解難題
1. 項(xiàng)目概述多峰高斯擬合的挑戰(zhàn)與價(jià)值在信號(hào)處理、光譜分析、色譜分離、生物醫(yī)學(xué)成像乃至金融數(shù)據(jù)分析中我們常常會(huì)遇到一種典型的數(shù)據(jù)形態(tài)一個(gè)看似復(fù)雜的波形實(shí)際上是由多個(gè)獨(dú)立的、相互重疊的“峰”疊加而成。比如一張質(zhì)譜圖上的多個(gè)離子峰一條光譜中不同元素的特征譜線(xiàn)或者心電圖里相鄰的P波、QRS波和T波。直接觀察這些混合在一起的峰我們很難精確地知道每個(gè)峰的中心位置、高度和寬度而這些參數(shù)恰恰是定量分析的核心。這時(shí)多峰高斯擬合就成了從混沌中提取秩序的“數(shù)學(xué)手術(shù)刀”。高斯函數(shù)或者說(shuō)正態(tài)分布曲線(xiàn)因其完美的鐘形對(duì)稱(chēng)性和良好的數(shù)學(xué)性質(zhì)成為描述這些獨(dú)立峰最常用的模型。單峰擬合相對(duì)簡(jiǎn)單但當(dāng)多個(gè)高斯峰擠在一起尤其是高度不同、寬度不一、甚至基線(xiàn)還有傾斜或偏移時(shí)問(wèn)題就變得棘手了。手動(dòng)“猜”參數(shù)幾乎不可能而簡(jiǎn)單的自動(dòng)擬合算法很容易陷入局部最優(yōu)解給出完全不合理的結(jié)果比如把兩個(gè)峰擬合成了一個(gè)寬峰或者擬合出的峰跑到數(shù)據(jù)范圍之外去了。我這次要分享的就是在MATLAB環(huán)境下成功實(shí)現(xiàn)三個(gè)重疊高斯峰的精確擬合全過(guò)程。這不僅僅是調(diào)用一個(gè)fit函數(shù)那么簡(jiǎn)單它涉及對(duì)數(shù)據(jù)本質(zhì)的理解、初始參數(shù)的巧妙估計(jì)、擬合算法的選擇以及大量“踩坑”后總結(jié)出的調(diào)試技巧。無(wú)論你是分析實(shí)驗(yàn)數(shù)據(jù)的研究生還是處理監(jiān)測(cè)數(shù)據(jù)的工程師這套從數(shù)據(jù)預(yù)處理、模型構(gòu)建、參數(shù)初始化到結(jié)果驗(yàn)證的完整流程都能讓你在面對(duì)復(fù)雜多峰數(shù)據(jù)時(shí)心里更有底。2. 核心思路與模型構(gòu)建理解“擬合”在做什么在動(dòng)手寫(xiě)代碼之前我們必須徹底想清楚多峰高斯擬合我們到底在求什么這決定了我們整個(gè)方案的架構(gòu)。2.1 數(shù)學(xué)模型三個(gè)高斯峰的疊加我們的目標(biāo)模型是三個(gè)高斯函數(shù)的線(xiàn)性疊加再加上一個(gè)可能存在的基線(xiàn)Baseline。一個(gè)標(biāo)準(zhǔn)的高斯函數(shù)公式如下y A * exp(-(x - μ)^2 / (2 * σ^2))其中A 峰高Amplitude。決定了峰的最大值。μ 峰位Mean/Center。決定了峰在x軸上的中心位置。σ 標(biāo)準(zhǔn)差Standard Deviation。決定了峰的寬度。半高全寬FWHM與σ的關(guān)系為FWHM 2√(2ln2) * σ ≈ 2.355 * σ。對(duì)于三個(gè)峰我們的總模型就是y_total Baseline Peak1 Peak2 Peak3即y_total (b0 b1*x) A1*exp(-(x-μ1)^2/(2*σ1^2)) A2*exp(-(x-μ2)^2/(2*σ2^2)) A3*exp(-(x-μ3)^2/(2*σ3^2))這里我特意將基線(xiàn)設(shè)為一次線(xiàn)性項(xiàng)b0 b1*x而不是一個(gè)常數(shù)。這是因?yàn)樵趯?shí)際數(shù)據(jù)中特別是光譜或色譜數(shù)據(jù)由于儀器背景或漂移基線(xiàn)傾斜非常常見(jiàn)。忽略它會(huì)導(dǎo)致峰高和峰面積的估計(jì)產(chǎn)生系統(tǒng)誤差。注意是否包含基線(xiàn)、基線(xiàn)是常數(shù)b0還是一次項(xiàng)b0b1*x甚至二次項(xiàng)需要根據(jù)你的數(shù)據(jù)實(shí)際情況判斷。一個(gè)簡(jiǎn)單的判斷方法是觀察數(shù)據(jù)中“無(wú)峰”區(qū)域的趨勢(shì)。如果拿不準(zhǔn)從簡(jiǎn)單模型常數(shù)基線(xiàn)開(kāi)始嘗試如果擬合殘差數(shù)據(jù)點(diǎn)與擬合曲線(xiàn)的差值呈現(xiàn)明顯的趨勢(shì)性分布則說(shuō)明需要更復(fù)雜的基線(xiàn)模型。2.2 擬合的本質(zhì)非線(xiàn)性最小二乘優(yōu)化擬合就是尋找一組模型參數(shù)上面提到的A1, μ1, σ1, A2, μ2, σ2, A3, μ3, σ3, b0, b1使得模型計(jì)算出的曲線(xiàn)y_total與實(shí)測(cè)數(shù)據(jù)點(diǎn)y_data之間的總體差異最小。這個(gè)差異通常用殘差平方和RSS來(lái)衡量RSS Σ(y_data_i - y_total_i)^2。擬合過(guò)程就是一個(gè)不斷調(diào)整這11個(gè)參數(shù)讓RSS達(dá)到最小的優(yōu)化過(guò)程。由于高斯函數(shù)是非線(xiàn)性的這是一個(gè)“非線(xiàn)性最小二乘”問(wèn)題。MATLAB的lsqcurvefit或曲線(xiàn)擬合工具箱的fit函數(shù)內(nèi)部使用的就是諸如“Levenberg-Marquardt”之類(lèi)的算法來(lái)解決這個(gè)問(wèn)題。這類(lèi)算法非常強(qiáng)大但有一個(gè)致命弱點(diǎn)高度依賴(lài)初始參數(shù)猜測(cè)。如果初始值離真實(shí)值太遠(yuǎn)算法極易收斂到錯(cuò)誤的局部最優(yōu)解而不是全局最優(yōu)解。因此整個(gè)多峰擬合成功的關(guān)鍵一半在于構(gòu)建正確的模型另一半就在于如何為這11個(gè)參數(shù)提供一個(gè)“聰明”的初始估計(jì)。接下來(lái)我們就進(jìn)入實(shí)戰(zhàn)環(huán)節(jié)。3. 數(shù)據(jù)準(zhǔn)備與初始參數(shù)估計(jì)為成功擬合奠基假設(shè)我們有一組實(shí)測(cè)數(shù)據(jù)x是自變量如波長(zhǎng)、時(shí)間、質(zhì)量數(shù)y是因變量如強(qiáng)度、吸光度、響應(yīng)值。數(shù)據(jù)已經(jīng)以數(shù)組形式存在于MATLAB工作區(qū)。3.1 數(shù)據(jù)可視化與初步觀察第一步永遠(yuǎn)是把數(shù)據(jù)畫(huà)出來(lái)用肉眼進(jìn)行第一次“診斷”。figure; plot(x, y, b.-, LineWidth, 1, MarkerSize, 10); xlabel(X (e.g., Wavelength)); ylabel(Y (e.g., Intensity)); title(Raw Data - Visual Inspection for Peaks); grid on;仔細(xì)觀察圖形識(shí)別峰的數(shù)量目標(biāo)是3個(gè)但你要確認(rèn)數(shù)據(jù)中是否明顯有3個(gè)凸起。有時(shí)噪聲或畸變會(huì)產(chǎn)生“假峰”。估計(jì)峰位μ用鼠標(biāo)光標(biāo)大致讀取三個(gè)峰頂對(duì)應(yīng)的x坐標(biāo)。記下它們例如mu1_guess, mu2_guess, mu3_guess。估計(jì)峰高A大致估計(jì)每個(gè)峰頂?shù)膟值。注意由于峰重疊這個(gè)值會(huì)低于該峰獨(dú)立存在時(shí)的高度?？梢韵扔梅逯禍p去附近“谷底”的y值來(lái)粗略估計(jì)。觀察基線(xiàn)看看數(shù)據(jù)最左邊和最右邊的點(diǎn)以及峰谷之間的區(qū)域連線(xiàn)趨勢(shì)是水平的還是傾斜的這決定了基線(xiàn)模型。3.2. 自動(dòng)化初始估計(jì)技巧對(duì)于更復(fù)雜或批量的數(shù)據(jù)我們可以借助MATLAB內(nèi)置函數(shù)進(jìn)行輔助估計(jì)這比肉眼估計(jì)更穩(wěn)健。1. 尋找峰值點(diǎn)估計(jì) μ 和 A使用findpeaks函數(shù)。這個(gè)函數(shù)能幫你找到局部極大值點(diǎn)并忽略一些小噪聲峰。[pks, locs, widths, proms] findpeaks(y, x, ... MinPeakProminence, max(y)*0.05, ... % 設(shè)置最小峰突出度過(guò)濾噪聲 MinPeakDistance, range(x)*0.05); % 設(shè)置最小峰間距避免識(shí)別到同一個(gè)峰pks是找到的峰值高度A的初始估計(jì)。locs是對(duì)應(yīng)的x位置μ的初始估計(jì)。widths和proms可以輔助估計(jì)σ。 檢查找到的峰數(shù)量是否大于等于3。如果太多可以調(diào)整MinPeakProminence和MinPeakDistance參數(shù)如果太少則調(diào)低這些閾值。2. 估計(jì)峰寬σ高斯峰的寬度σ可以通過(guò)多種方式估計(jì)半高寬法對(duì)于一個(gè)孤立的峰找到峰高一半處的兩個(gè)點(diǎn)其x坐標(biāo)差值即為FWHM然后除以2.355得到σ。對(duì)于重疊峰這很難自動(dòng)完成。利用findpeaks的輸出findpeaks函數(shù)返回的widths是每個(gè)峰在半高處的寬度這近似就是FWHM。我們可以用sigma_guess widths / 2.355。經(jīng)驗(yàn)公式如果數(shù)據(jù)點(diǎn)間隔均勻可以觀察峰上升沿/下降沿的陡峭程度。一個(gè)粗略的起點(diǎn)是設(shè)σ的初始值為相鄰峰間距的1/5到1/10。3. 估計(jì)基線(xiàn)參數(shù)b0, b1一個(gè)簡(jiǎn)單有效的方法是對(duì)數(shù)據(jù)兩端例如前5%和后5%的數(shù)據(jù)點(diǎn)進(jìn)行線(xiàn)性擬合得到的截距和斜率作為b0和b1的初始值。n length(x); indices [1:round(n*0.05), round(n*0.95):n]; % 取頭尾5%的索引 p polyfit(x(indices), y(indices), 1); % 一階多項(xiàng)式擬合 b0_guess p(2); b1_guess p(1);假設(shè)通過(guò)以上方法我們得到了如下初始猜測(cè)% 峰1 (最左側(cè)) A1_guess pks(1); mu1_guess locs(1); sigma1_guess widths(1)/2.355; % 峰2 (中間) A2_guess pks(2); mu2_guess locs(2); sigma2_guess widths(2)/2.355; % 峰3 (最右側(cè)) A3_guess pks(3); mu3_guess locs(3); sigma3_guess widths(3)/2.355; % 基線(xiàn) % b0_guess, b1_guess 來(lái)自polyfit將這些初始值組合成一個(gè)向量initial_guess [A1_guess, mu1_guess, sigma1_guess, A2_guess, mu2_guess, sigma2_guess, A3_guess, mu3_guess, sigma3_guess, b0_guess, b1_guess]。4. 擬合實(shí)現(xiàn)兩種主流方法與詳細(xì)步驟有了模型和初始參數(shù)我們就可以開(kāi)始擬合了。這里介紹兩種最常用的方法使用lsqcurvefit優(yōu)化函數(shù)和使用曲線(xiàn)擬合工具箱的fit函數(shù)。lsqcurvefit更底層、靈活適合集成到腳本中fit函數(shù)更直觀、快捷適合交互式分析。4.1 方法一使用lsqcurvefit進(jìn)行擬合lsqcurvefit是優(yōu)化工具箱中的函數(shù)它直接處理非線(xiàn)性最小二乘問(wèn)題。第一步定義模型函數(shù)在MATLAB中創(chuàng)建一個(gè)函數(shù)文件例如multiGauss.m或者使用匿名函數(shù)。這里用匿名函數(shù)示例% 定義三峰高斯帶線(xiàn)性基線(xiàn)的模型函數(shù) % params: [A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1] multiGaussModel (params, x) ... params(1) * exp(-(x - params(2)).^2 / (2 * params(3)^2)) ... % 峰1 params(4) * exp(-(x - params(41)).^2 / (2 * params(42)^2)) ... % 峰2 (注意索引) params(7) * exp(-(x - params(71)).^2 / (2 * params(72)^2)) ... % 峰3 params(10) params(11) * x; % 線(xiàn)性基線(xiàn) b0 b1*x注意索引的對(duì)應(yīng)關(guān)系。為了清晰也可以將參數(shù)解包multiGaussModel (p, x) ... p(1)*exp(-(x-p(2)).^2/(2*p(3)^2)) ... p(4)*exp(-(x-p(5)).^2/(2*p(6)^2)) ... p(7)*exp(-(x-p(8)).^2/(2*p(9)^2)) ... p(10) p(11)*x;第二步設(shè)置邊界約束關(guān)鍵步驟這是避免擬合出荒謬結(jié)果如負(fù)的峰寬、峰位跑到天涯海角的關(guān)鍵。我們需要為每個(gè)參數(shù)設(shè)置合理的上下界lb和ub。% 基于初始猜測(cè)設(shè)置邊界 % 順序: [A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1] % 下界 (Lower Bounds) lb [0, min(x), 0, ... % 峰1: 振幅0, 峰位在數(shù)據(jù)范圍內(nèi) 寬度0 0, min(x), 0, ... % 峰2 0, min(x), 0, ... % 峰3 -inf, -inf]; % 基線(xiàn)參數(shù)可以為任意值 % 上界 (Upper Bounds) ub [inf, max(x), range(x)/2, ... % 峰1: 振幅無(wú)上限峰位在數(shù)據(jù)范圍內(nèi)寬度小于數(shù)據(jù)范圍一半合理假設(shè) inf, max(x), range(x)/2, ... % 峰2 inf, max(x), range(x)/2, ... % 峰3 inf, inf]; % 基線(xiàn)參數(shù)無(wú)限制 % 可以更精細(xì)地約束例如讓峰位按順序排列避免擬合時(shí)峰位互換 % lb(2) lb(5) lb(8) 且 ub(2) ub(5) ub(8) 的邏輯更復(fù)雜通?？亢玫某跏贾当苊?。第三步執(zhí)行擬合% 設(shè)置優(yōu)化選項(xiàng)提高顯示細(xì)節(jié) options optimoptions(lsqcurvefit, Display, iter, Algorithm, trust-region-reflective); % ‘levenberg-marquardt’算法不支持邊界這里用‘trust-region-reflective’ % 執(zhí)行擬合 [params_fitted, resnorm, residual, exitflag, output] ... lsqcurvefit(multiGaussModel, initial_guess, x, y, lb, ub, options); disp(擬合參數(shù) (A, mu, sigma, b0, b1):); disp(params_fitted);params_fitted就是擬合得到的最優(yōu)參數(shù)向量。4.2 方法二使用曲線(xiàn)擬合工具箱fit函數(shù)fit函數(shù)語(yǔ)法更貼近“擬合”這個(gè)概念并且能自動(dòng)生成豐富的統(tǒng)計(jì)信息和繪圖。第一步定義擬合類(lèi)型和選項(xiàng)% 使用 fittype 定義模型coefficients 指定參數(shù)名稱(chēng) ft fittype(A1*exp(-(x-mu1)^2/(2*sigma1^2)) A2*exp(-(x-mu2)^2/(2*sigma2^2)) A3*exp(-(x-mu3)^2/(2*sigma3^2)) b0 b1*x, ... independent, x, ... dependent, y, ... coefficients, {A1, mu1, sigma1, A2, mu2, sigma2, A3, mu3, sigma3, b0, b1}); % 設(shè)置擬合選項(xiàng)包括初始值和邊界 opts fitoptions(ft); opts.StartPoint initial_guess; % 傳入我們之前準(zhǔn)備好的初始猜測(cè)向量 opts.Lower lb; % 下界 opts.Upper ub; % 上界 opts.Display Iter; % 顯示迭代過(guò)程 % opts.Robust LAR; % 如果數(shù)據(jù)有異常點(diǎn)可以嘗試穩(wěn)健擬合第二步執(zhí)行擬合并繪圖% 執(zhí)行擬合 [fitresult, gof] fit(x, y, ft, opts); % 顯示擬合結(jié)果 disp(fitresult); disp(gof); % 輸出擬合優(yōu)度統(tǒng)計(jì)量如 R-square, RMSE % 繪制擬合結(jié)果 figure; plot(fitresult, x, y); legend(原始數(shù)據(jù), 擬合曲線(xiàn), Location, Best); xlabel(X); ylabel(Y); title(三峰高斯擬合結(jié)果);fitresult是一個(gè)包含所有擬合參數(shù)的對(duì)象可以通過(guò)fitresult.A1,fitresult.mu1等方式訪問(wèn)。gof包含了擬合優(yōu)度的信息如決定系數(shù)rsquare越接近1說(shuō)明擬合越好。5. 結(jié)果評(píng)估、可視化與問(wèn)題排查擬合完成并不意味著結(jié)束我們必須嚴(yán)格評(píng)估擬合質(zhì)量并診斷可能的問(wèn)題。5.1 可視化評(píng)估四象限診斷圖一張好的診斷圖勝過(guò)千言萬(wàn)語(yǔ)。我習(xí)慣同時(shí)繪制四個(gè)子圖figure(Position, [100, 100, 1200, 800]); % 子圖1原始數(shù)據(jù) vs. 擬合曲線(xiàn) subplot(2,2,1); plot(x, y, b., MarkerSize, 8); hold on; x_fine linspace(min(x), max(x), 1000); % 生成更密的點(diǎn)用于繪制光滑曲線(xiàn) y_fitted multiGaussModel(params_fitted, x_fine); % 或用 fitresult(x_fine) plot(x_fine, y_fitted, r-, LineWidth, 2); legend(Data, Fitted Curve, Location, Best); title(Fit Overview); grid on; % 子圖2殘差圖 (Residuals) subplot(2,2,2); residuals y - multiGaussModel(params_fitted, x); % 計(jì)算殘差 plot(x, residuals, k^, MarkerSize, 5, MarkerFaceColor, k); hold on; plot([min(x), max(x)], [0,0], r--); % 零參考線(xiàn) xlabel(X); ylabel(Residual); title(Residual Plot); grid on; % 好的擬合殘差應(yīng)隨機(jī)分布在零線(xiàn)附近無(wú)趨勢(shì)性。 % 子圖3殘差直方圖 subplot(2,2,3); histogram(residuals, 20, Normalization, probability, FaceColor, [0.5, 0.5, 0.5]); xlabel(Residual); ylabel(Probability); title(Residual Distribution); grid on; % 理想情況應(yīng)接近均值為0的正態(tài)分布。 % 子圖4分峰顯示 (Peak Decomposition) subplot(2,2,4); plot(x_fine, y_fitted, k-, LineWidth, 1.5); hold on; % 計(jì)算并繪制每個(gè)單獨(dú)的峰 peak1 params_fitted(1) * exp(-(x_fine - params_fitted(2)).^2 / (2 * params_fitted(3)^2)); peak2 params_fitted(4) * exp(-(x_fine - params_fitted(5)).^2 / (2 * params_fitted(6)^2)); peak3 params_fitted(7) * exp(-(x_fine - params_fitted(8)).^2 / (2 * params_fitted(9)^2)); baseline params_fitted(10) params_fitted(11) * x_fine; plot(x_fine, peak1, g--, LineWidth, 1); plot(x_fine, peak2, b--, LineWidth, 1); plot(x_fine, peak3, m--, LineWidth, 1); plot(x_fine, baseline, c:, LineWidth, 1); legend(Total Fit, Peak 1, Peak 2, Peak 3, Baseline, Location, Best); title(Peak Decomposition); grid on;通過(guò)這四張圖你可以一目了然地判斷總覽圖擬合曲線(xiàn)是否完美貼合數(shù)據(jù)點(diǎn)殘差圖殘差是否隨機(jī)、無(wú)規(guī)律如果呈現(xiàn)“U”型或“∩”型說(shuō)明模型選擇不當(dāng)如基線(xiàn)模型不對(duì)。殘差分布是否近似正態(tài)嚴(yán)重偏離可能暗示有異常點(diǎn)或模型系統(tǒng)誤差。分峰圖分解出的單個(gè)峰是否合理有沒(méi)有出現(xiàn)負(fù)峰圖形上表現(xiàn)為向下凸峰位是否與預(yù)期相符5.2 定量評(píng)估指標(biāo)除了看圖還要看數(shù)決定系數(shù) R2gof.rsquare。大于0.99通常說(shuō)明擬合很好但要注意對(duì)于非常尖銳的峰即使R2很高峰面積也可能不準(zhǔn)。殘差平方和 (RSS) 或 均方根誤差 (RMSE)sqrt(gof.sse / length(x))。越小越好但要在不同數(shù)據(jù)集間比較才有意義。參數(shù)置信區(qū)間使用confint(fitresult)可以計(jì)算參數(shù)的95%置信區(qū)間。如果某個(gè)參數(shù)的置信區(qū)間非常寬例如包含0或負(fù)值說(shuō)明該參數(shù)不可靠可能是數(shù)據(jù)信息不足或者該參數(shù)與其他參數(shù)強(qiáng)相關(guān)共線(xiàn)性。5.3 常見(jiàn)問(wèn)題與排查技巧實(shí)錄在實(shí)際操作中你幾乎一定會(huì)遇到下面這些問(wèn)題。以下是我的排查清單問(wèn)題1擬合失敗提示“未收斂”或“達(dá)到最大迭代次數(shù)”。原因初始值太差或者邊界設(shè)置不合理導(dǎo)致優(yōu)化算法找不到下降方向。解決放松邊界先將所有邊界設(shè)得非常寬如lb -inf(1,11); ub inf(1,11);只保留sigma0這樣的物理約束。如果能擬合再逐步收緊邊界。改進(jìn)初始值回到第3步用更穩(wěn)健的方法如對(duì)數(shù)據(jù)平滑后再找峰估計(jì)初始值?？梢試L試手動(dòng)在圖上選點(diǎn)。分步擬合先擬合一個(gè)峰固定其參數(shù)再加入第二個(gè)峰擬合以此類(lèi)推。這能有效降低優(yōu)化難度。換用算法lsqcurvefit可以嘗試‘levenberg-marquardt’算法但不支持邊界。fit函數(shù)可以嘗試‘Robust’選項(xiàng)。問(wèn)題2擬合結(jié)果中某個(gè)峰的振幅是負(fù)值或者峰寬極大/極小。原因典型的局部最優(yōu)解或者模型過(guò)于復(fù)雜過(guò)擬合。解決施加物理約束強(qiáng)制振幅A大于0峰寬σ在一個(gè)合理范圍內(nèi)如數(shù)據(jù)范圍的1/100到1/2。簡(jiǎn)化模型檢查是否真的需要三個(gè)峰也許兩個(gè)峰加一個(gè)更復(fù)雜的基線(xiàn)模型就夠了?；蛘邍L試固定其中一兩個(gè)你認(rèn)為最確定的參數(shù)如已知某個(gè)峰位。檢查數(shù)據(jù)質(zhì)量數(shù)據(jù)噪聲是否太大考慮先對(duì)原始數(shù)據(jù)進(jìn)行平滑處理如Savitzky-Golay濾波但注意平滑可能扭曲峰形。問(wèn)題3殘差圖顯示出明顯的規(guī)律性如彎曲趨勢(shì)。原因模型不足以描述數(shù)據(jù)。通常是基線(xiàn)模型不合適。解決升級(jí)基線(xiàn)模型將常數(shù)基線(xiàn)b0改為線(xiàn)性b0b1*x甚至二次b0b1*xb2*x^2。檢查峰函數(shù)模型數(shù)據(jù)峰形是否不對(duì)稱(chēng)高斯函數(shù)是對(duì)稱(chēng)的。如果峰有明顯拖尾可能需要考慮洛倫茲Lorentzian函數(shù)或二者的混合Voigt profile。此時(shí)模型應(yīng)改為A / (1 ((x-μ)/σ)^2)。問(wèn)題4兩個(gè)峰的參數(shù)擬合后幾乎一樣或者峰位互換了。原因初始值中兩個(gè)峰的估計(jì)位置太接近或者算法在迭代中發(fā)生了“跳變”。解決嚴(yán)格約束峰位順序在lsqcurvefit中可以通過(guò)設(shè)置非線(xiàn)性約束來(lái)實(shí)現(xiàn)但這比較復(fù)雜。更實(shí)用的方法是在初始值中明確指定mu1_guess mu2_guess mu3_guess并設(shè)置不重疊的邊界如[mu1_guess-Δ, mu1_guessΔ]。使用“鎖定”策略先擬合最左側(cè)和最右側(cè)的兩個(gè)峰將它們的參數(shù)固定再擬合中間的峰。問(wèn)題5擬合速度很慢尤其是數(shù)據(jù)點(diǎn)很多的時(shí)候。原因每次迭代都要計(jì)算整個(gè)模型函數(shù)數(shù)據(jù)點(diǎn)多則計(jì)算量大。解決數(shù)據(jù)降采樣在保持峰形特征的前提下對(duì)數(shù)據(jù)進(jìn)行適當(dāng)降采樣。提供解析雅可比矩陣對(duì)于lsqcurvefit可以提供一個(gè)函數(shù)來(lái)計(jì)算模型關(guān)于各個(gè)參數(shù)的導(dǎo)數(shù)雅可比矩陣這能極大加速收斂。但對(duì)于高斯模型手動(dòng)推導(dǎo)并編寫(xiě)雅可比矩陣比較繁瑣除非對(duì)性能有極致要求否則通常不需要。6. 進(jìn)階技巧與擴(kuò)展應(yīng)用當(dāng)你掌握了三峰擬合后可以嘗試以下進(jìn)階操作讓分析更上一層樓。6.1 自動(dòng)化與批處理如果你有成百上千條光譜需要分析手動(dòng)操作是不可行的。你需要將上述流程封裝成函數(shù)。function [fittedParams, gofStats, fitResult] fitThreeGaussPeaks(xData, yData, initialGuess, lowerBounds, upperBounds) % 封裝三峰高斯擬合流程 % 輸入xData, yData, 初始猜測(cè)下界上界 % 輸出擬合參數(shù)統(tǒng)計(jì)量fit結(jié)果對(duì)象 ft fittype(...); % 定義模型 opts fitoptions(...); % 設(shè)置選項(xiàng) % ... [填充具體代碼] [fitResult, gofStats] fit(xData, yData, ft, opts); fittedParams coeffvalues(fitResult); end然后在一個(gè)循環(huán)中調(diào)用這個(gè)函數(shù)處理每個(gè)數(shù)據(jù)文件。關(guān)鍵難點(diǎn)在于自動(dòng)生成可靠的初始猜測(cè)。你可以開(kāi)發(fā)一個(gè)穩(wěn)健的峰值檢測(cè)和基線(xiàn)估計(jì)子函數(shù)作為fitThreeGaussPeaks的前置步驟。6.2 從擬合參數(shù)到物理量峰面積計(jì)算在很多應(yīng)用中峰高A受儀器條件影響大而峰面積Area更能代表物質(zhì)的量。對(duì)于高斯峰面積S A * σ * √(2π)。A1 params_fitted(1); sigma1 params_fitted(3); area1 A1 * sigma1 * sqrt(2*pi); % 同理計(jì)算 area2, area3如果存在線(xiàn)性基線(xiàn)上述公式計(jì)算的是峰相對(duì)于基線(xiàn)的凈面積這是正確的。6.3 模型選擇高斯 vs. 洛倫茲 vs. 沃伊特不是所有的峰都是完美的高斯形。高斯峰源于多普勒增寬、某些色譜過(guò)程。峰形較“瘦”衰減更快。洛倫茲峰源于自然增寬、某些共振現(xiàn)象。峰形較“胖”有更長(zhǎng)的拖尾。沃伊特峰高斯和洛倫茲的卷積能描述更復(fù)雜的增寬機(jī)制。在MATLAB中只需修改fittype中的表達(dá)式即可切換模型洛倫茲‘A / (1 ((x-mu)/sigma)^2)’沃伊特需要自定義函數(shù)或使用Faddeeva函數(shù)近似較為復(fù)雜。選擇模型的依據(jù)一是物理過(guò)程的先驗(yàn)知識(shí)二是看哪種模型的殘差更小、更隨機(jī)。可以都試一下用gof.rsquare和殘差圖來(lái)輔助判斷。6.4 不確定性分析與誤差傳遞擬合出的參數(shù)是有不確定性的置信區(qū)間。當(dāng)我們用這些參數(shù)計(jì)算衍生量如峰面積、峰位差時(shí)誤差也會(huì)傳遞。 MATLAB的曲線(xiàn)擬合工具箱可以計(jì)算參數(shù)的雅可比矩陣和協(xié)方差矩陣。對(duì)于簡(jiǎn)單的線(xiàn)性函數(shù)如面積計(jì)算誤差傳遞可以用公式近似。但對(duì)于復(fù)雜情況推薦使用蒙特卡洛模擬假設(shè)擬合參數(shù)服從以最佳估計(jì)值為均值、以標(biāo)準(zhǔn)誤差為方差的多維正態(tài)分布。從這個(gè)分布中隨機(jī)抽取大量如10000組參數(shù)集。對(duì)每組參數(shù)計(jì)算你關(guān)心的衍生量如面積。衍生量結(jié)果的分布如直方圖就給出了該量的估計(jì)值及其不確定性如95%置信區(qū)間。這個(gè)過(guò)程在MATLAB中實(shí)現(xiàn)起來(lái)需要一些編程但它提供了最全面的不確定性評(píng)估。成功擬合三個(gè)重疊高斯峰是一個(gè)從理論到實(shí)踐再到經(jīng)驗(yàn)積累的完整過(guò)程。它考驗(yàn)的不僅僅是對(duì)MATLAB函數(shù)的熟悉程度更是對(duì)數(shù)據(jù)、模型和優(yōu)化算法的綜合理解。最深刻的體會(huì)是沒(méi)有“一鍵萬(wàn)能”的擬合。初始值的精心設(shè)置、物理約束的合理施加、以及基于殘差圖的模型診斷這些手動(dòng)步驟的重要性往往超過(guò)選擇哪個(gè)具體的擬合函數(shù)。每次擬合都是一次與數(shù)據(jù)的對(duì)話(huà)你需要不斷提出問(wèn)題模型對(duì)嗎初始值好嗎并根據(jù)數(shù)據(jù)的“回答”殘差圖、參數(shù)值來(lái)調(diào)整策略。當(dāng)你看到分解出的三個(gè)光滑鐘形曲線(xiàn)完美地拼合成原始數(shù)據(jù)時(shí)那種從雜亂中提煉出清晰信息的成就感正是數(shù)據(jù)分析工作最大的樂(lè)趣所在。