實(shí)戰(zhàn):從擬合到異常檢測)
1. 為什么搞數(shù)據(jù)統(tǒng)計(jì)繞不開概率分布做數(shù)據(jù)分析這些年我發(fā)現(xiàn)一個挺有意思的現(xiàn)象很多人拿到一組數(shù)據(jù)第一反應(yīng)就是算均值、算方差然后畫個直方圖看看形狀覺得這就是統(tǒng)計(jì)分析了。但真正做過幾個完整項(xiàng)目之后你會發(fā)現(xiàn)均值方差只是最表層的描述真正決定你分析結(jié)論靠不靠譜的是你對數(shù)據(jù)背后概率分布的理解。舉個我踩過的坑。早些年做一個傳感器信號處理的項(xiàng)目采集了幾萬條數(shù)據(jù)我直接用正態(tài)分布假設(shè)去算置信區(qū)間結(jié)果發(fā)現(xiàn)尾部數(shù)據(jù)完全對不上后來一查才知道那批數(shù)據(jù)是典型的重尾分布用正態(tài)分布去套置信區(qū)間嚴(yán)重偏窄導(dǎo)致誤判率飆升。從那以后我就養(yǎng)成了一個習(xí)慣拿到任何數(shù)據(jù)先看它服從什么分布再談后續(xù)分析。這篇內(nèi)容主要聊兩件事一是Matlab里常見的概率分布怎么用、怎么選二是分位數(shù)這個被很多人低估但極其好用的工具。適合已經(jīng)會用Matlab做基礎(chǔ)數(shù)據(jù)處理、想往統(tǒng)計(jì)分析方向深入的朋友也適合做信號處理、金融分析、質(zhì)量控制、機(jī)器學(xué)習(xí)特征工程的同學(xué)參考。不需要你有多深的數(shù)學(xué)功底但基本的概率論概念得知道比如什么是概率密度、什么是累積分布。我會把每個分布的核心參數(shù)、Matlab里的調(diào)用方式、實(shí)際場景怎么選、以及分位數(shù)在工程中的妙用都講清楚盡量做到你讀完就能直接上手改代碼。2. 常見概率分布的核心邏輯與選型思路2.1 分布選型的底層邏輯先看數(shù)據(jù)怎么產(chǎn)生的很多人學(xué)概率分布是死記硬背公式我覺得這是最沒效率的方式。正確的思路是先搞清楚數(shù)據(jù)是怎么產(chǎn)生的分布自然就出來了。打個比方正態(tài)分布為什么這么常見因?yàn)樗嵌鄠€獨(dú)立小因素疊加的結(jié)果。比如測量誤差是溫度波動、儀器精度、讀數(shù)偏差等一堆小因素加在一起造成的根據(jù)中心極限定理這種疊加就趨向正態(tài)。而指數(shù)分布描述的是等待時間比如一個設(shè)備從開始運(yùn)行到出故障的時間這種無記憶性的等待過程就是指數(shù)分布。所以選分布的時候我會問自己三個問題這個數(shù)據(jù)是計(jì)數(shù)還是連續(xù)測量計(jì)數(shù)通常往泊松、二項(xiàng)靠連續(xù)測量往正態(tài)、指數(shù)、伽馬靠。數(shù)據(jù)有沒有上下界比如比例數(shù)據(jù)在0到1之間那貝塔分布可能比正態(tài)更合適。數(shù)據(jù)的尾部厚不厚金融收益率、網(wǎng)絡(luò)流量這類數(shù)據(jù)尾部往往很厚正態(tài)會低估極端事件。這三個問題想清楚了分布選型基本不會跑偏。2.2 Matlab里分布對象的統(tǒng)一調(diào)用范式Matlab從R2013a開始引入了概率分布對象這套體系我覺得這是它做統(tǒng)計(jì)分析最舒服的地方。所有分布都遵循同一套命名規(guī)則學(xué)會一個就會用全部。核心函數(shù)就四個pdf/cdf概率密度/累積分布icdf逆累積分布也就是給定概率求分位數(shù)random生成服從該分布的隨機(jī)數(shù)fit/mle從數(shù)據(jù)擬合分布參數(shù)命名規(guī)則是分布名pdf、分布名cdf這種形式。比如正態(tài)是normpdf、normcdf、norminv指數(shù)是exppdf、expcdf、expinv。注意逆函數(shù)有的叫inv有的叫icdf新版統(tǒng)一推薦用對象方式調(diào)用。我一般這樣用% 創(chuàng)建正態(tài)分布對象 pd makedist(Normal, mu, 0, sigma, 1); % 求x1.96處的累積概率 p cdf(pd, 1.96); % 求95%分位數(shù) q icdf(pd, 0.95); % 生成1000個隨機(jī)數(shù) r random(pd, 1000, 1);這種對象式寫法的好處是參數(shù)只寫一次后面所有操作都復(fù)用代碼干凈很多。以前用normcdf(x, mu, sigma)這種函數(shù)式寫法每個調(diào)用都要重復(fù)傳參數(shù)參數(shù)一多就容易出錯。2.3 幾個高頻分布的實(shí)戰(zhàn)場景對照我把工程中最常用的幾個分布整理成一張表方便你快速對照分布核心參數(shù)典型場景Matlab對象名正態(tài)分布mu, sigma測量誤差、身高體重、信號噪聲Normal指數(shù)分布mu設(shè)備壽命、等待時間、到達(dá)間隔Exponential泊松分布lambda單位時間事件數(shù)、缺陷數(shù)Poisson二項(xiàng)分布n, p合格率、點(diǎn)擊率、成功次數(shù)Binomial伽馬分布a, b累積等待時間、降雨量Gamma貝塔分布a, b比例、概率的先驗(yàn)Beta威布爾分布a, b可靠性壽命、風(fēng)速Weibull對數(shù)正態(tài)分布mu, sigma收入、股價、粒徑Lognormal這張表我建議你存下來做項(xiàng)目的時候?qū)χ?。特別提醒一點(diǎn)指數(shù)分布是伽馬分布的特例形狀參數(shù)為1威布爾分布又比指數(shù)更靈活能描述遞增或遞減的故障率。做可靠性分析的時候威布爾往往比指數(shù)更貼合實(shí)際因?yàn)楹芏嘣O(shè)備的故障率不是恒定的。2.4 參數(shù)估計(jì)別只會用均值方差選好分布之后下一步是從數(shù)據(jù)里估計(jì)參數(shù)。這里有個坑很多人直接用樣本均值和樣本方差去套正態(tài)分布參數(shù)這在數(shù)據(jù)量大的時候沒問題但數(shù)據(jù)量小或者有離群點(diǎn)的時候會偏得厲害。Matlab提供了fitdist函數(shù)用最大似然估計(jì)來擬合參數(shù)比矩估計(jì)更穩(wěn)健% 假設(shè)data是你的樣本數(shù)據(jù) pd fitdist(data, Normal); % 查看估計(jì)出的參數(shù) pd.mu pd.sigma % 擬合威布爾分布 pdw fitdist(data, Weibull);fitdist還支持指定分布類型后自動選最優(yōu)參數(shù)甚至可以用fitdist配合AIC準(zhǔn)則在多個候選分布里挑最合適的。我一般會同時擬合正態(tài)、威布爾、對數(shù)正態(tài)三個然后比較AIC值選最小的那個。這個流程在可靠性分析里幾乎是標(biāo)配。注意fitdist對數(shù)據(jù)的預(yù)處理有要求如果數(shù)據(jù)里有NaN或者Inf會直接報錯。擬合前記得用data data(isfinite(data))清理一下。3. 分位數(shù)被低估的統(tǒng)計(jì)分析利器3.1 分位數(shù)到底是什么為什么比均值好用分位數(shù)這個概念說白了就是把數(shù)據(jù)從小到大排好隊(duì)然后按比例切幾刀。中位數(shù)就是50%分位數(shù)把數(shù)據(jù)一分為二四分位數(shù)把數(shù)據(jù)切成四段百分位數(shù)切成一百段。為什么我說它比均值好用因?yàn)榫祵﹄x群點(diǎn)極其敏感而分位數(shù)幾乎不受影響。舉個例子一個部門10個人9個人月薪1萬1個人月薪100萬平均下來每人10.9萬但這個均值完全不能代表這個部門的真實(shí)收入水平。而中位數(shù)還是1萬穩(wěn)穩(wěn)地反映了大多數(shù)人的情況。在工程上分位數(shù)的應(yīng)用場景特別多質(zhì)量控制用P9999%分位數(shù)來設(shè)定告警閾值比用均值加幾倍標(biāo)準(zhǔn)差更靠譜因?yàn)轫憫?yīng)時間這類數(shù)據(jù)往往不是正態(tài)的。金融風(fēng)控VaR在險價值本質(zhì)上就是某個置信水平下的分位數(shù)。性能分析看接口響應(yīng)時間的P50、P95、P99比看平均值有意義得多。異常檢測用IQR四分位距來判定離群點(diǎn)這是箱線圖的底層邏輯。3.2 Matlab里算分位數(shù)的幾種方式Matlab算分位數(shù)主要有三個函數(shù)各有適用場景第一種quantile函數(shù)最直接data randn(1000, 1); % 求0.25, 0.5, 0.75分位數(shù) q quantile(data, [0.25, 0.5, 0.75]); % 求P99 p99 quantile(data, 0.99);quantile默認(rèn)用的是線性插值法對中間分位數(shù)比較友好。但它有個細(xì)節(jié)要注意當(dāng)分位數(shù)落在兩個數(shù)據(jù)點(diǎn)之間時不同的插值方法結(jié)果不一樣。Matlab支持多種方法通過Method參數(shù)指定比如exact、approximate、midpoint等。數(shù)據(jù)量小的時候建議用exact數(shù)據(jù)量大的時候用默認(rèn)的近似方法更快。第二種prctile函數(shù)專門算百分位數(shù)% 求第90百分位數(shù) p90 prctile(data, 90);這個函數(shù)和quantile本質(zhì)一樣只是輸入用百分?jǐn)?shù)而不是小數(shù)看個人習(xí)慣。第三種分布對象的icdf當(dāng)你已知分布類型時用pd fitdist(data, Normal); % 求理論上的95%分位數(shù) q95 icdf(pd, 0.95);這種方式的好處是即使你的樣本里沒有極端值也能外推出理論分位數(shù)。比如你只有1000個樣本想估計(jì)P99.9直接算樣本分位數(shù)會很不準(zhǔn)但用擬合的分布去算就靠譜得多。3.3 分位數(shù)與分布的關(guān)系一張圖看懂分位數(shù)和累積分布函數(shù)CDF其實(shí)是一體兩面的關(guān)系。CDF是給定x求概率分位數(shù)是給定概率求x兩者互為逆函數(shù)。我經(jīng)常用這個關(guān)系來做可視化。比如畫一條CDF曲線然后在縱軸上標(biāo)出0.25、0.5、0.75對應(yīng)的橫軸值就是四分位數(shù)。這樣一眼就能看出數(shù)據(jù)的分布形態(tài)如果CDF曲線在中間很陡說明數(shù)據(jù)集中如果很平說明分散。data exprnd(2, 1000, 1); [f, x] ecdf(data); plot(x, f, LineWidth, 1.5); hold on; % 標(biāo)出中位數(shù) med quantile(data, 0.5); plot([med med], [0 0.5], r--); plot([0 med], [0.5 0.5], r--); xlabel(數(shù)值); ylabel(累積概率); title(經(jīng)驗(yàn)CDF與中位數(shù)); grid on;這段代碼畫出來的圖能直觀展示分位數(shù)在分布中的位置。做報告的時候放這么一張圖比干巴巴的數(shù)字有說服力得多。3.4 分位數(shù)在異常檢測中的實(shí)戰(zhàn)用法這是我最想強(qiáng)調(diào)的一個應(yīng)用。很多人做異常檢測還在用均值±3倍標(biāo)準(zhǔn)差這套但這個方法有個致命前提數(shù)據(jù)必須近似正態(tài)。一旦數(shù)據(jù)偏態(tài)或者有重尾這個方法就會失效。用IQR法則就穩(wěn)健得多Q1 quantile(data, 0.25); Q3 quantile(data, 0.75); IQR Q3 - Q1; lower_bound Q1 - 1.5 * IQR; upper_bound Q3 1.5 * IQR; outliers data(data lower_bound | data upper_bound);這個1.5倍IQR的系數(shù)是經(jīng)驗(yàn)值來自箱線圖的定義。它的好處是不依賴分布假設(shè)對偏態(tài)數(shù)據(jù)也能用。如果數(shù)據(jù)特別重尾可以把系數(shù)調(diào)到3減少誤報。我在一個日志分析項(xiàng)目里就用過這套。當(dāng)時要檢測異常的請求延遲用均值加標(biāo)準(zhǔn)差的方法正常波動都被判成異常了換成IQR法則之后誤報率直接降了一個數(shù)量級。這個經(jīng)驗(yàn)我后來在好幾個項(xiàng)目里復(fù)用屢試不爽。4. 完整實(shí)操從數(shù)據(jù)到分布擬合再到分位數(shù)分析4.1 數(shù)據(jù)準(zhǔn)備與探索性分析假設(shè)我們拿到一批設(shè)備運(yùn)行時長數(shù)據(jù)存在lifetime變量里。第一步永遠(yuǎn)是探索性分析別急著套模型。% 清理無效值 lifetime lifetime(isfinite(lifetime) lifetime 0); % 基礎(chǔ)統(tǒng)計(jì)量 fprintf(樣本量: %d\n, numel(lifetime)); fprintf(均值: %.2f\n, mean(lifetime)); fprintf(中位數(shù): %.2f\n, median(lifetime)); fprintf(標(biāo)準(zhǔn)差: %.2f\n, std(lifetime)); fprintf(偏度: %.2f\n, skewness(lifetime)); fprintf(峰度: %.2f\n, kurtosis(lifetime));這里重點(diǎn)看偏度和峰度。偏度接近0說明對稱大于0說明右偏長尾在右邊小于0說明左偏。峰度接近3說明接近正態(tài)大于3說明尖峰重尾。設(shè)備壽命數(shù)據(jù)通常是右偏的偏度為正這時候用正態(tài)分布就不合適了。再看直方圖和Q-Q圖subplot(1,2,1); histogram(lifetime, 50, Normalization, pdf); hold on; subplot(1,2,2); qqplot(lifetime);Q-Q圖是判斷分布類型的神器。如果數(shù)據(jù)點(diǎn)大致落在一條直線上說明服從對應(yīng)的分布。Matlab的qqplot默認(rèn)對比正態(tài)分布如果點(diǎn)明顯彎曲就說明不是正態(tài)。4.2 多分布擬合與AIC準(zhǔn)則選優(yōu)探索完之后我會同時擬合幾個候選分布然后用AIC挑最優(yōu)dist_names {Normal, Exponential, Weibull, Lognormal, Gamma}; aic_values zeros(size(dist_names)); for i 1:numel(dist_names) try pd fitdist(lifetime, dist_names{i}); aic_values(i) aic(pd); fprintf(%s: AIC %.2f\n, dist_names{i}, aic_values(i)); catch aic_values(i) Inf; fprintf(%s: 擬合失敗\n, dist_names{i}); end end [~, best_idx] min(aic_values); fprintf(最優(yōu)分布: %s\n, dist_names{best_idx});AIC的規(guī)則是越小越好它同時考慮了擬合優(yōu)度和模型復(fù)雜度避免過擬合。這個流程我建議你封裝成一個函數(shù)以后每個項(xiàng)目直接調(diào)用。實(shí)操心得擬合失敗最常見的原因是數(shù)據(jù)里有零或負(fù)值而指數(shù)、威布爾、伽馬、對數(shù)正態(tài)都要求數(shù)據(jù)為正。擬合前一定要檢查數(shù)據(jù)范圍必要時做平移或者換分布。4.3 參數(shù)解讀與擬合優(yōu)度檢驗(yàn)選好分布之后要解讀參數(shù)的實(shí)際含義。以威布爾分布為例它有兩個參數(shù)形狀參數(shù)a和尺度參數(shù)b。a 1故障率遞減說明早期故障多產(chǎn)品在磨合期。a 1故障率恒定退化成指數(shù)分布。a 1故障率遞增說明產(chǎn)品在老化期。這個解讀在可靠性工程里非常有用能直接指導(dǎo)維護(hù)策略。如果a 1說明設(shè)備越用越容易壞該考慮定期更換了。擬合優(yōu)度檢驗(yàn)用卡方檢驗(yàn)或者K-S檢驗(yàn)pd fitdist(lifetime, Weibull); % K-S檢驗(yàn) [h, p] kstest(lifetime, CDF, pd); if h 0 fprintf(不能拒絕原假設(shè)擬合可接受p %.4f\n, p); else fprintf(擬合不佳p %.4f\n, p); endK-S檢驗(yàn)的原假設(shè)是數(shù)據(jù)服從該分布p值大于0.05就不能拒絕說明擬合可接受。但要注意樣本量很大的時候K-S檢驗(yàn)會過于敏感一點(diǎn)點(diǎn)偏差都會判為不擬合。這時候要結(jié)合Q-Q圖綜合判斷別只看p值。4.4 基于擬合分布計(jì)算關(guān)鍵分位數(shù)擬合好分布之后就可以算各種分位數(shù)了。這里我一般會同時算經(jīng)驗(yàn)分位數(shù)和理論分位數(shù)做對比pd fitdist(lifetime, Weibull); probs [0.01, 0.05, 0.25, 0.5, 0.75, 0.95, 0.99]; emp_q quantile(lifetime, probs); theo_q icdf(pd, probs); T table(probs, emp_q, theo_q, ... VariableNames, {Probability, Empirical, Theoretical}); disp(T);如果兩者差距不大說明擬合得不錯如果尾部差距大說明分布選得不對得回去重新選。在實(shí)際項(xiàng)目里我最常用的是P95和P99。比如設(shè)定設(shè)備維護(hù)閾值用P95意味著95%的情況下不會誤報只有真正異常的5%才會觸發(fā)。這個閾值比用均值加標(biāo)準(zhǔn)差靈活得多因?yàn)槟憧梢愿鶕?jù)業(yè)務(wù)容忍度調(diào)整分位水平。5. 常見問題與排查技巧實(shí)錄5.1 分布擬合報錯怎么辦這是新手最常遇到的問題。我把常見報錯和解決方法整理成表報錯信息原因解決方法Data must be positive數(shù)據(jù)含零或負(fù)值檢查數(shù)據(jù)剔除或平移NaN or Inf detected數(shù)據(jù)含無效值data data(isfinite(data))Not enough data樣本量太小至少30個樣本最好100以上Parameter estimate failed數(shù)據(jù)退化全相同檢查數(shù)據(jù)是否常量我踩過最坑的一次是數(shù)據(jù)里混了幾個負(fù)值因?yàn)閭鞲衅髋紶枙祷劐e誤碼。擬合一直失敗查了半天才發(fā)現(xiàn)。后來我養(yǎng)成了習(xí)慣擬合前先跑一遍summary(data)把數(shù)據(jù)的min、max、NaN數(shù)量都看一眼能省很多調(diào)試時間。5.2 分位數(shù)計(jì)算結(jié)果和預(yù)期不符有時候你會發(fā)現(xiàn)quantile算出來的結(jié)果和手算的不一樣這通常是插值方法的問題。Matlab的quantile默認(rèn)對中間分位數(shù)用線性插值對極端分位數(shù)用其他方法。如果你想要和某些統(tǒng)計(jì)軟件比如R的type 7一致的結(jié)果需要顯式指定方法q quantile(data, 0.5, Method, exact);不同方法在數(shù)據(jù)量大的時候差異很小但數(shù)據(jù)量小的時候可能差挺多。做嚴(yán)謹(jǐn)分析的時候建議在報告里注明用的哪種方法。5.3 理論分位數(shù)和經(jīng)驗(yàn)分位數(shù)差很多這個問題的根源通常是分布選錯了。比如數(shù)據(jù)明明是重尾的你硬套正態(tài)分布那P99的理論值會明顯小于經(jīng)驗(yàn)值。排查思路是先畫Q-Q圖看尾部是否偏離再對比幾個候選分布的AIC最后用K-S檢驗(yàn)確認(rèn)。三步下來基本能定位問題。還有一種可能是樣本量太小經(jīng)驗(yàn)分位數(shù)本身就不穩(wěn)定。比如只有50個樣本算P99基本沒有意義因?yàn)?9%分位對應(yīng)的數(shù)據(jù)點(diǎn)根本不在樣本里。這時候只能依賴?yán)碚摲治粩?shù)但要清楚它的不確定性很大。5.4 實(shí)操避坑清單最后分享幾條我總結(jié)的避坑經(jīng)驗(yàn)都是真金白銀換來的永遠(yuǎn)先畫圖再建模。直方圖、Q-Q圖、箱線圖三件套能幫你避開80%的坑。分位數(shù)比均值穩(wěn)健但不是萬能。數(shù)據(jù)量小于30的時候分位數(shù)波動很大別太當(dāng)真。擬合分布要交叉驗(yàn)證。把數(shù)據(jù)分成訓(xùn)練集和測試集在測試集上驗(yàn)證擬合效果避免過擬合。AIC不是唯一標(biāo)準(zhǔn)。業(yè)務(wù)可解釋性有時候比AIC更重要選一個業(yè)務(wù)上說得通的分布比選一個AIC最小但沒人懂的分布強(qiáng)。記錄你的分析過程。分布選型、參數(shù)、分位數(shù)方法都記下來下次復(fù)用或者別人接手的時候能省大量時間。這套流程我在多個項(xiàng)目里反復(fù)打磨過從傳感器數(shù)據(jù)分析到用戶行為分析都用得上。核心就一句話先理解數(shù)據(jù)怎么產(chǎn)生的再選分布最后用分位數(shù)做決策。順序別搞反了搞反了就是本末倒置。如果你剛開始接觸這塊建議拿一組自己熟悉的數(shù)據(jù)練手把整個流程走一遍比看十篇文章都管用。