1. 項(xiàng)目概述與核心問(wèn)題拆解去年帶隊(duì)參加高教社杯數(shù)學(xué)建模競(jìng)賽C題“古代玻璃制品的成分分析與鑒別”給我留下了深刻印象。這道題之所以經(jīng)典是因?yàn)樗昝赖貙⒊橄蟮臄?shù)學(xué)模型與具體的文化遺產(chǎn)保護(hù)問(wèn)題結(jié)合了起來(lái)考察的遠(yuǎn)不止是編程和算法更是對(duì)數(shù)據(jù)本質(zhì)的理解和解決實(shí)際問(wèn)題的邏輯構(gòu)建能力。題目提供了一批古代玻璃文物的化學(xué)成分檢測(cè)數(shù)據(jù)要求我們通過(guò)這些數(shù)據(jù)判斷文物的類型、風(fēng)化情況并分析其成分間的關(guān)聯(lián)規(guī)律。這聽起來(lái)像是一個(gè)標(biāo)準(zhǔn)的分類和回歸問(wèn)題但當(dāng)你真正深入數(shù)據(jù)時(shí)會(huì)發(fā)現(xiàn)處處是坑數(shù)據(jù)有缺失、成分間存在嚴(yán)重的多重共線性、風(fēng)化過(guò)程并非簡(jiǎn)單的線性變化。如果只是機(jī)械地套用幾個(gè)機(jī)器學(xué)習(xí)模型大概率會(huì)得到一個(gè)看似漂亮但實(shí)際毫無(wú)解釋力的結(jié)果。這篇分享我就結(jié)合我們當(dāng)時(shí)的解題思路和代碼實(shí)現(xiàn)來(lái)詳細(xì)拆解如何系統(tǒng)性地應(yīng)對(duì)這類成分分析問(wèn)題希望能給正在備戰(zhàn)數(shù)模或?qū)?shù)據(jù)分析感興趣的朋友一些實(shí)在的參考。簡(jiǎn)單來(lái)說(shuō)這個(gè)題目的核心目標(biāo)有三個(gè)第一根據(jù)化學(xué)成分對(duì)玻璃文物進(jìn)行類型鑒別比如是高鉀玻璃還是鉛鋇玻璃第二分析風(fēng)化過(guò)程判斷文物表面是否風(fēng)化并研究風(fēng)化對(duì)成分的影響第三探索成分之間的關(guān)聯(lián)規(guī)律為文物研究和保護(hù)提供依據(jù)。數(shù)據(jù)是典型的“寬表”即樣本數(shù)不多但每個(gè)樣本的特征化學(xué)成分氧化物含量很多且這些特征之和理論上應(yīng)為100%即“定和約束”這直接限制了許多統(tǒng)計(jì)方法的應(yīng)用。因此整個(gè)解題過(guò)程就是一場(chǎng)與數(shù)據(jù)特性斗智斗勇的旅程。2. 解題核心思路與整體方案設(shè)計(jì)面對(duì)這樣的問(wèn)題一個(gè)清晰的、分階段的解決方案至關(guān)重要。我們的整體思路可以概括為“數(shù)據(jù)預(yù)處理奠基統(tǒng)計(jì)探索先行模型選擇求精結(jié)果解釋為魂”。很多新手團(tuán)隊(duì)容易犯的錯(cuò)誤是一上來(lái)就想著用最復(fù)雜的模型比如深度學(xué)習(xí)去擬合忽略了數(shù)據(jù)本身的特質(zhì)和問(wèn)題的物理背景最終導(dǎo)致模型過(guò)擬合或結(jié)果無(wú)法解釋。2.1 第一階段數(shù)據(jù)預(yù)處理與探索性數(shù)據(jù)分析這是所有數(shù)據(jù)分析項(xiàng)目的基石對(duì)于成分?jǐn)?shù)據(jù)尤其關(guān)鍵。原始數(shù)據(jù)通常存在缺失值、異常值并且成分?jǐn)?shù)據(jù)是“閉合數(shù)據(jù)”。2.1.1 缺失值處理策略數(shù)據(jù)中部分化學(xué)成分存在缺失。對(duì)于成分?jǐn)?shù)據(jù)不能簡(jiǎn)單地用均值或中位數(shù)填充因?yàn)檫@會(huì)破壞數(shù)據(jù)的“定和”結(jié)構(gòu)所有成分百分比之和應(yīng)為100%。我們采用了兩種策略對(duì)于個(gè)別零星缺失如果某個(gè)成分在大多數(shù)樣本中均被檢測(cè)到僅在個(gè)別樣本中缺失且該成分含量通常較低我們考慮使用多重插補(bǔ)方法在保持成分?jǐn)?shù)據(jù)協(xié)方差結(jié)構(gòu)的前提下進(jìn)行填充。在Python中可以使用IterativeImputer。對(duì)于大量缺失或未檢測(cè)成分如果某個(gè)成分如某些微量元素在大量樣本中均為“未檢測(cè)”我們將其視為一種“結(jié)構(gòu)性零值”。在后續(xù)分析中我們有兩種選擇一是將該成分從特征集中移除避免引入過(guò)多噪聲二是在進(jìn)行對(duì)數(shù)比變換時(shí)采用一個(gè)極小的值如檢測(cè)限的一半進(jìn)行替換但這需要謹(jǐn)慎因?yàn)闀?huì)影響變換的穩(wěn)定性。2.1.2 成分?jǐn)?shù)據(jù)的特殊性處理——對(duì)數(shù)比變換這是本題最核心的統(tǒng)計(jì)知識(shí)點(diǎn)之一。普通的歐氏距離、相關(guān)性系數(shù)對(duì)于成分?jǐn)?shù)據(jù)是失效的因?yàn)槌煞种g不是獨(dú)立的一個(gè)成分的增加必然導(dǎo)致其他成分的減少。我們必須將數(shù)據(jù)從“單純形空間”映射到“歐氏空間”。最常用的方法是中心對(duì)數(shù)比變換。import numpy as np import pandas as pd def clr_transform(data): 對(duì)成分?jǐn)?shù)據(jù)進(jìn)行中心對(duì)數(shù)比變換。 data: DataFrame, 每一行是一個(gè)樣本每一列是一種成分比例或百分比。 # 確保數(shù)據(jù)為正將0替換為一個(gè)極小值需根據(jù)實(shí)際情況調(diào)整 data data.replace(0, 1e-6) # 計(jì)算幾何均值 gmean np.exp(np.mean(np.log(data), axis1)) # CLR變換 clr_data np.log(data.div(gmean, axis0)) return clr_data經(jīng)過(guò)CLR變換后的數(shù)據(jù)各個(gè)特征近似服從正態(tài)分布且消除了定和約束此時(shí)才能安全地使用PCA、聚類、回歸等基于歐氏距離的多元統(tǒng)計(jì)方法。2.1.3 異常值檢測(cè)我們使用馬氏距離結(jié)合主成分分析后的得分來(lái)檢測(cè)多元異常值。在CLR變換后的數(shù)據(jù)上計(jì)算馬氏距離距離過(guò)大的樣本可能測(cè)量有誤或?qū)儆谔厥忸悇e需要結(jié)合文物背景知識(shí)進(jìn)行甄別決定是保留、修正還是剔除。注意數(shù)據(jù)預(yù)處理的所有步驟都必須記錄在案并在論文中詳細(xì)說(shuō)明理由。評(píng)委非常看重處理過(guò)程的科學(xué)性和嚴(yán)謹(jǐn)性。2.2 第二階段統(tǒng)計(jì)分析與可視化探索在正式建模前通過(guò)統(tǒng)計(jì)圖表直觀理解數(shù)據(jù)能為我們后續(xù)的模型選擇提供至關(guān)重要的方向。描述性統(tǒng)計(jì)與箱線圖計(jì)算各類玻璃高鉀、鉛鋇各成分含量的均值、中位數(shù)、標(biāo)準(zhǔn)差、極值。繪制箱線圖可以直觀對(duì)比兩類玻璃在主要成分如K2O, PbO, BaO上的分布差異初步判斷哪些成分是關(guān)鍵的鑒別指標(biāo)。相關(guān)性分析熱圖計(jì)算CLR變換后數(shù)據(jù)中各成分之間的皮爾遜相關(guān)系數(shù)并用熱圖可視化。這能幫助我們理解成分間的共生或拮抗關(guān)系。例如PbO和BaO可能呈現(xiàn)強(qiáng)正相關(guān)這與鉛鋇玻璃的工藝特點(diǎn)是吻合的。主成分分析散點(diǎn)圖對(duì)CLR變換后的數(shù)據(jù)做PCA將前兩個(gè)或三個(gè)主成分的得分繪制成散點(diǎn)圖并用文物類型或風(fēng)化狀態(tài)著色。如果散點(diǎn)圖上不同類型能明顯分開說(shuō)明成分?jǐn)?shù)據(jù)本身具有很好的區(qū)分度簡(jiǎn)單的線性分類器如LDA就可能有效。3. 核心模型構(gòu)建與代碼實(shí)現(xiàn)詳解基于探索性分析的結(jié)果我們針對(duì)三個(gè)子問(wèn)題構(gòu)建了相應(yīng)的模型。3.1 子問(wèn)題一玻璃類型鑒別分類問(wèn)題這是一個(gè)典型的二分類問(wèn)題高鉀 vs. 鉛鋇。我們采用了模型集成的思路不依賴單一模型以提升穩(wěn)健性。3.1.1 特征選擇并非所有化學(xué)成分都是有效的鑒別特征。我們使用遞歸特征消除結(jié)合交叉驗(yàn)證來(lái)選擇最具判別力的成分。from sklearn.feature_selection import RFECV from sklearn.linear_model import LogisticRegression from sklearn.model_selection import StratifiedKFold # X_train_clr 是CLR變換后的訓(xùn)練特征 y_train 是類型標(biāo)簽 lr LogisticRegression(max_iter1000, solverliblinear) rfecv RFECV(estimatorlr, step1, cvStratifiedKFold(5), scoringaccuracy) rfecv.fit(X_train_clr, y_train) print(最優(yōu)特征數(shù)量, rfecv.n_features_) print(被選中的特征, X_train_clr.columns[rfecv.support_])3.1.2 分類模型構(gòu)建與比較我們訓(xùn)練了多個(gè)分類器并在獨(dú)立的驗(yàn)證集上比較性能邏輯回歸可解釋性強(qiáng)能給出特征系數(shù)判斷哪些成分對(duì)分類貢獻(xiàn)大。支持向量機(jī)適用于小樣本、高維數(shù)據(jù)特別是當(dāng)兩類數(shù)據(jù)在特征空間里非線性可分時(shí)高斯核SVM往往有不錯(cuò)效果。隨機(jī)森林能自動(dòng)處理特征交互給出特征重要性排序?qū)Ξ惓Ｖ挡幻舾小GBoost/LightGBM梯度提升樹通常能獲得最高的準(zhǔn)確率但需要仔細(xì)調(diào)參以防過(guò)擬合。from sklearn.svm import SVC from sklearn.ensemble import RandomForestClassifier from sklearn.model_selection import cross_val_score models { LR: LogisticRegression(max_iter1000), SVM: SVC(kernelrbf, probabilityTrue), RF: RandomForestClassifier(n_estimators100, random_state42) } for name, model in models.items(): scores cross_val_score(model, X_train_selected, y_train, cv5, scoringaccuracy) print(f{name} 平均交叉驗(yàn)證準(zhǔn)確率{scores.mean():.4f} (/- {scores.std()*2:.4f}))3.1.3 模型集成與最終預(yù)測(cè)我們選擇表現(xiàn)最好的2-3個(gè)模型采用軟投票法進(jìn)行集成。即每個(gè)模型對(duì)測(cè)試樣本預(yù)測(cè)其屬于各類的概率然后將這些概率平均取概率最高的類別作為最終預(yù)測(cè)結(jié)果。這種方法通常比硬投票直接投票類別或單一模型更穩(wěn)定。from sklearn.ensemble import VotingClassifier voting_clf VotingClassifier( estimators[(lr, best_lr), (svm, best_svm), (rf, best_rf)], votingsoft # 軟投票 ) voting_clf.fit(X_train_selected, y_train) final_predictions voting_clf.predict(X_test_selected)3.2 子問(wèn)題二風(fēng)化情況分析與影響研究這部分是問(wèn)題的難點(diǎn)和亮點(diǎn)。我們不能僅僅做一個(gè)“是否風(fēng)化”的分類更要定量研究風(fēng)化過(guò)程。3.2.1 風(fēng)化程度量化題目只給了“風(fēng)化”與“無(wú)風(fēng)化”的標(biāo)簽但實(shí)際風(fēng)化是一個(gè)連續(xù)過(guò)程。我們創(chuàng)新性地嘗試構(gòu)建一個(gè)風(fēng)化指數(shù)。思路是風(fēng)化會(huì)導(dǎo)致某些易溶成分如K2O, Na2O流失而某些穩(wěn)定成分如SiO2, Al2O3相對(duì)富集。我們可以通過(guò)主成分分析找到最能區(qū)分風(fēng)化與未風(fēng)化樣本的方向即第一主成分并將樣本在該主成分上的得分標(biāo)準(zhǔn)化后作為“風(fēng)化指數(shù)”。指數(shù)越高代表風(fēng)化特征越顯著。3.2.2 風(fēng)化影響的統(tǒng)計(jì)分析差異性檢驗(yàn)對(duì)每個(gè)化學(xué)成分使用Mann-Whitney U檢驗(yàn)非參數(shù)檢驗(yàn)對(duì)分布假設(shè)要求低比較其在風(fēng)化與未風(fēng)化組間的中位數(shù)是否存在顯著差異。找出風(fēng)化過(guò)程中顯著流失或富集的成分。相關(guān)性分析計(jì)算各成分含量與上述“風(fēng)化指數(shù)”之間的斯皮爾曼秩相關(guān)系數(shù)進(jìn)一步驗(yàn)證其與風(fēng)化程度的相關(guān)性。可視化繪制堆疊面積圖或平行坐標(biāo)圖直觀展示從無(wú)風(fēng)化到嚴(yán)重風(fēng)化按風(fēng)化指數(shù)排序的樣本其成分比例的變化趨勢(shì)。3.2.3 風(fēng)化前后成分變化的預(yù)測(cè)模型這是一個(gè)回歸問(wèn)題根據(jù)未風(fēng)化玻璃的成分預(yù)測(cè)其風(fēng)化后的表面成分。我們采用了以下步驟數(shù)據(jù)配對(duì)理想情況下需要有同一器物風(fēng)化前后對(duì)應(yīng)的數(shù)據(jù)。本題數(shù)據(jù)可能不具備嚴(yán)格配對(duì)我們可以假設(shè)同類型、成分相近的未風(fēng)化玻璃其風(fēng)化行為相似。因此可以為每個(gè)風(fēng)化樣本在未風(fēng)化樣本中尋找其K近鄰基于CLR變換后的成分?jǐn)?shù)據(jù)用這些近鄰的未風(fēng)化成分均值作為其“風(fēng)化前”的估計(jì)。建立回歸模型以估計(jì)的“風(fēng)化前”成分為自變量以實(shí)際測(cè)得的“風(fēng)化后”表面成分為因變量為每個(gè)受風(fēng)化影響的成分如K2O建立回歸模型如嶺回歸、彈性網(wǎng)絡(luò)以處理多重共線性。解讀模型回歸系數(shù)可以解釋為在其他成分不變的情況下該成分風(fēng)化前后的保留比例或變化規(guī)律。3.3 子問(wèn)題三成分關(guān)聯(lián)與亞類劃分這部分旨在探索數(shù)據(jù)內(nèi)部更深層的結(jié)構(gòu)。3.3.1 相關(guān)性網(wǎng)絡(luò)分析基于CLR變換后數(shù)據(jù)的相關(guān)系數(shù)矩陣我們可以構(gòu)建一個(gè)相關(guān)性網(wǎng)絡(luò)。將每種成分視為一個(gè)節(jié)點(diǎn)如果兩種成分之間的相關(guān)系數(shù)絕對(duì)值超過(guò)一個(gè)閾值如0.7就在它們之間連一條邊。然后使用社區(qū)發(fā)現(xiàn)算法如Louvain算法對(duì)網(wǎng)絡(luò)進(jìn)行劃分得到的社區(qū)可能就是具有共同地球化學(xué)行為或工藝來(lái)源的成分組合。這比單純看熱圖更直觀。3.3.2 無(wú)監(jiān)督聚類發(fā)現(xiàn)亞類在高鉀玻璃和鉛鋇玻璃內(nèi)部是否還存在不同的配方或產(chǎn)地亞類我們使用聚類算法來(lái)探索。方法對(duì)每一大類玻璃的CLR數(shù)據(jù)分別進(jìn)行譜聚類或DBSCAN聚類。譜聚類能發(fā)現(xiàn)非球形的簇DBSCAN能自動(dòng)識(shí)別噪聲點(diǎn)。驗(yàn)證聚類結(jié)果需要結(jié)合輪廓系數(shù)和實(shí)際意義來(lái)判斷。例如如果聚類出的亞類在PbO/BaO比值、或微量元素組合上表現(xiàn)出明顯差異這很可能對(duì)應(yīng)不同的制作時(shí)期或作坊。可視化使用t-SNE或UMAP將高維數(shù)據(jù)降至2維或3維進(jìn)行可視化并用聚類標(biāo)簽著色可以非常清晰地展示亞類的存在。4. 代碼實(shí)現(xiàn)中的關(guān)鍵技巧與避坑指南紙上談兵終覺(jué)淺代碼實(shí)現(xiàn)時(shí)才會(huì)遇到真正的挑戰(zhàn)。這里分享幾個(gè)我們踩過(guò)坑后總結(jié)的關(guān)鍵點(diǎn)。4.1 成分?jǐn)?shù)據(jù)中零值的處理這是最大的坑之一。CLR變換需要取對(duì)數(shù)零值會(huì)導(dǎo)致無(wú)窮大。簡(jiǎn)單地用一個(gè)固定小值如1e-6替換所有零值可能會(huì)嚴(yán)重扭曲低含量成分的信息特別是當(dāng)數(shù)據(jù)中存在大量“未檢測(cè)”時(shí)。我們的經(jīng)驗(yàn)是首先區(qū)分“真零”工藝上就不含該成分和“檢測(cè)限下的零”實(shí)際有但儀器未檢出。這需要領(lǐng)域知識(shí)。對(duì)于“檢測(cè)限下的零”可以用多重插補(bǔ)法或者使用專門為成分?jǐn)?shù)據(jù)設(shè)計(jì)的Coda方法中的zCompositions包R語(yǔ)言或scikit-composition庫(kù)Python早期版本進(jìn)行處理。在純Python環(huán)境中一種穩(wěn)健的做法是采用貝葉斯多重插補(bǔ)假設(shè)數(shù)據(jù)服從狄利克雷分布。4.2 避免“數(shù)據(jù)泄露”在構(gòu)建風(fēng)化預(yù)測(cè)模型時(shí)尋找未風(fēng)化樣本的K近鄰必須在訓(xùn)練集內(nèi)進(jìn)行絕對(duì)不能用測(cè)試集的數(shù)據(jù)來(lái)為訓(xùn)練集樣本找近鄰。整個(gè)流程應(yīng)該是劃分訓(xùn)練集和測(cè)試集按樣本劃分。在訓(xùn)練集內(nèi)部為每個(gè)風(fēng)化樣本找其K近鄰來(lái)自訓(xùn)練集中的未風(fēng)化樣本構(gòu)建回歸模型。用訓(xùn)練好的模型去預(yù)測(cè)測(cè)試集中的風(fēng)化樣本。 這個(gè)順序一旦搞反就會(huì)得到過(guò)于樂(lè)觀的、不可信的結(jié)果。4.3 模型的可解釋性優(yōu)先數(shù)學(xué)建模競(jìng)賽不是單純的機(jī)器學(xué)習(xí)比賽模型結(jié)果需要能說(shuō)得通。例如在分類模型中如果邏輯回歸顯示PbO的系數(shù)極大且為正那么我們可以合理地解釋為“PbO含量高是鉛鋇玻璃的決定性特征”。但如果一個(gè)復(fù)雜的深度學(xué)習(xí)模型達(dá)到了同樣的準(zhǔn)確率卻無(wú)法提供這樣的解釋其得分可能反而不如邏輯回歸。因此在模型選擇上我們遵循“簡(jiǎn)單有效可解釋”優(yōu)先的原則。4.4 結(jié)果的穩(wěn)健性檢驗(yàn)對(duì)于關(guān)鍵結(jié)論一定要做穩(wěn)健性檢驗(yàn)。例如更換隨機(jī)種子在涉及隨機(jī)性的步驟如數(shù)據(jù)劃分、隨機(jī)森林、K-Means聚類中多次運(yùn)行觀察結(jié)果是否穩(wěn)定。擾動(dòng)輸入數(shù)據(jù)對(duì)成分?jǐn)?shù)據(jù)加入微小的高斯噪聲看分類準(zhǔn)確率或聚類結(jié)果是否發(fā)生劇烈變化。使用不同的預(yù)處理方法比如對(duì)比CLR變換和等距對(duì)數(shù)比變換的結(jié)果差異。5. 論文寫作與結(jié)果呈現(xiàn)要點(diǎn)模型和代碼只是基礎(chǔ)如何將你的工作清晰、有說(shuō)服力地呈現(xiàn)出來(lái)才是贏得評(píng)委青睞的關(guān)鍵。5.1 圖文并茂一圖勝千言使用組合圖例如將PCA散點(diǎn)圖、成分箱線圖和相關(guān)性熱圖組合在一張圖上分面展示不同類型、不同風(fēng)化狀態(tài)下的數(shù)據(jù)特征。流程圖必不可少繪制一張清晰的算法或分析流程圖讓評(píng)委一眼看懂你的技術(shù)路線。結(jié)果對(duì)比可視化對(duì)于分類結(jié)果一定要畫混淆矩陣對(duì)于聚類結(jié)果用輪廓圖展示聚類質(zhì)量。5.2 量化表述避免模糊不要說(shuō)“模型效果很好”要說(shuō)“集成模型在測(cè)試集上的準(zhǔn)確率達(dá)到94.2%F1-score為0.93”。不要說(shuō)“A成分和B成分有關(guān)”要說(shuō)“A與B的斯皮爾曼相關(guān)系數(shù)為0.82p0.01呈強(qiáng)正相關(guān)”。在分析風(fēng)化影響時(shí)給出具體的流失百分比估算如“風(fēng)化導(dǎo)致K2O平均流失約60%”。5.3 討論部分體現(xiàn)深度這是區(qū)分優(yōu)秀論文和普通論文的地方。不能只羅列結(jié)果要解釋結(jié)果背后的原因。聯(lián)系化學(xué)與考古學(xué)知識(shí)例如解釋為什么K2O在風(fēng)化中易流失鉀離子易溶于水為什么PbO和BaO高度相關(guān)因?yàn)樗鼈兂Ｒ糟U鋇礦的形式一起加入作為助熔劑。分析模型的局限性誠(chéng)實(shí)地指出你的方法基于哪些假設(shè)如風(fēng)化過(guò)程均勻、數(shù)據(jù)缺失是隨機(jī)的這些假設(shè)如果被違反會(huì)對(duì)結(jié)論產(chǎn)生什么影響。提出進(jìn)一步研究建議基于你的發(fā)現(xiàn)可以提出哪些新的問(wèn)題例如“本文發(fā)現(xiàn)的鉛鋇玻璃亞類一和亞類二是否對(duì)應(yīng)于不同的歷史時(shí)期建議結(jié)合出土墓葬年代信息進(jìn)行驗(yàn)證。”回顧整個(gè)解題過(guò)程最大的體會(huì)是數(shù)學(xué)建模競(jìng)賽考察的是一種系統(tǒng)性的問(wèn)題解決能力。從理解一個(gè)跨領(lǐng)域考古化學(xué)的實(shí)際問(wèn)題到將其轉(zhuǎn)化為數(shù)學(xué)語(yǔ)言再到選擇合適的工具進(jìn)行處理最后將數(shù)學(xué)結(jié)果翻譯回實(shí)際結(jié)論每一步都需要嚴(yán)謹(jǐn)?shù)乃伎己筒粩嗟臋?quán)衡。代碼能力很重要但它只是實(shí)現(xiàn)想法的工具。更核心的是對(duì)數(shù)據(jù)特性的洞察、對(duì)模型假設(shè)的理解以及將復(fù)雜結(jié)果清晰傳達(dá)的能力。這道玻璃成分分析題就是一個(gè)絕佳的練兵場(chǎng)它教會(huì)我們面對(duì)真實(shí)世界紛繁復(fù)雜的數(shù)據(jù)沒(méi)有銀彈模型只有具體問(wèn)題具體分析的審慎和創(chuàng)造力。