1. 項(xiàng)目背景與問(wèn)題引入當(dāng)數(shù)學(xué)模型遇上公共衛(wèi)生危機(jī)2022年初我參與了一場(chǎng)特殊的數(shù)學(xué)建模競(jìng)賽——認(rèn)證杯SPSSPRO杯。當(dāng)時(shí)全球仍籠罩在新冠疫情的陰影下我們團(tuán)隊(duì)抽到的C題第二階段直接將我們拉入了一個(gè)極具現(xiàn)實(shí)意義的戰(zhàn)場(chǎng)如何利用“污水流行病學(xué)”的原理為疫情防控提供量化的決策支持這個(gè)題目讓我印象深刻因?yàn)樗昝赖卦忈屃藬?shù)學(xué)建模如何從象牙塔走向街頭巷尾解決真實(shí)世界的棘手問(wèn)題。你可能聽(tīng)說(shuō)過(guò)通過(guò)檢測(cè)污水來(lái)追蹤毒品濫用情況這就是污水流行病學(xué)的經(jīng)典應(yīng)用。其核心邏輯很簡(jiǎn)單人體代謝后的許多物質(zhì)會(huì)隨排泄物進(jìn)入城市污水系統(tǒng)。通過(guò)檢測(cè)污水中這些特定生物標(biāo)志物比如病毒RNA片段、藥物殘留的濃度結(jié)合污水流量、服務(wù)人口等數(shù)據(jù)就能反推出一個(gè)社區(qū)內(nèi)某種物質(zhì)的消費(fèi)量或某種病原體的流行情況。這種方法具有匿名、客觀、覆蓋面廣、可早期預(yù)警的巨大優(yōu)勢(shì)。在新冠疫情中它被用來(lái)監(jiān)測(cè)社區(qū)感染水平甚至早于臨床報(bào)告發(fā)現(xiàn)疫情苗頭。然而題目給我們的不是一篇綜述而是一系列冰冷的數(shù)據(jù)和幾個(gè)尖銳的問(wèn)題。它模擬了一個(gè)城市的污水監(jiān)測(cè)數(shù)據(jù)要求我們建立數(shù)學(xué)模型估算真實(shí)感染人數(shù)分析數(shù)據(jù)中的異常波動(dòng)并評(píng)估不同防控措施的效果。這聽(tīng)起來(lái)像是一個(gè)完美的“黑箱”預(yù)測(cè)問(wèn)題但實(shí)際操作中從污水濃度到感染人數(shù)中間隔著人口流動(dòng)、病毒載量差異、檢測(cè)靈敏度、管道沉降、水力停留時(shí)間等無(wú)數(shù)個(gè)“灰箱”。題目提供的“校正因子”關(guān)鍵詞正是打開(kāi)這些灰箱的一把關(guān)鍵鑰匙。整個(gè)解題過(guò)程實(shí)際上是一場(chǎng)關(guān)于“如何用數(shù)學(xué)語(yǔ)言描述并量化不確定性”的思維體操。2. 核心問(wèn)題拆解從污水?dāng)?shù)據(jù)到管理決策的三重挑戰(zhàn)面對(duì)題目給出的數(shù)據(jù)集和問(wèn)題描述我們首先需要摒棄“直接套公式”的幻想。污水流行病學(xué)應(yīng)用于新冠監(jiān)測(cè)其建模鏈條比傳統(tǒng)應(yīng)用更為復(fù)雜。我們將賽題的核心訴求拆解為三個(gè)環(huán)環(huán)相扣的層次這構(gòu)成了我們整個(gè)論文的骨架。2.1 第一層反演問(wèn)題——如何從濃度估算感染人數(shù)這是最根本的數(shù)學(xué)模型問(wèn)題。輸入是污水中新冠病毒RNA的濃度單位拷貝數(shù)/升輸出是監(jiān)測(cè)區(qū)域內(nèi)可能的活躍感染人數(shù)。這個(gè)過(guò)程絕非簡(jiǎn)單的線性比例計(jì)算。我們需要構(gòu)建一個(gè)“反演模型”?；竟娇蚣苋缦耓 N_{infected} \frac{C_{wastewater} \times F \times V}{f \times R \times S \times \eta} ]其中(N_{infected})估算的感染人數(shù)。(C_{wastewater})污水中測(cè)得的病毒RNA濃度拷貝數(shù)/L。(F)污水稀釋因子考慮地下水滲入、雨水等。(V)污水廠日均處理流量L/天。(f)每人每天平均糞便排放量L/天/人。這是一個(gè)生理常數(shù)但個(gè)體差異大。(R)感染者日均通過(guò)糞便排出的病毒RNA拷貝數(shù)。這是最不確定、最核心的參數(shù)它取決于病程階段、個(gè)體免疫狀況、病毒毒株等多種因素。(S)病毒RNA在污水輸送和處理過(guò)程中的存活率/回收率。病毒RNA會(huì)降解也會(huì)被吸附在管道污泥上。(\eta)實(shí)驗(yàn)室檢測(cè)方法的回收率/效率。題目中提到的“校正因子”本質(zhì)上就是針對(duì) (R)、(S)、(\eta) 等這些高度不確定的參數(shù)進(jìn)行整體性或分項(xiàng)式的標(biāo)定與修正。我們的首要任務(wù)就是設(shè)計(jì)一種算法或模型來(lái)合理地確定或校準(zhǔn)這個(gè)綜合校正因子 (K 1/(f \times R \times S \times \eta))使得模型估算結(jié)果與已知的、有限的其他數(shù)據(jù)如同期臨床報(bào)告病例數(shù)在統(tǒng)計(jì)意義上最為吻合。注意這里切忌直接使用文獻(xiàn)中的單一數(shù)值。不同城市、不同季節(jié)、不同毒株、不同污水管網(wǎng)材質(zhì)這些參數(shù)都可能不同。建模的價(jià)值就在于根據(jù)本地?cái)?shù)據(jù)動(dòng)態(tài)地確定或校準(zhǔn)這些參數(shù)。2.2 第二層時(shí)空分析——如何識(shí)別和解釋數(shù)據(jù)異常題目提供的污水監(jiān)測(cè)數(shù)據(jù)通常包含時(shí)間序列每日/每周濃度可能還有空間序列來(lái)自城市不同采樣點(diǎn)的濃度。第二層挑戰(zhàn)在于分析這些數(shù)據(jù)的波動(dòng)模式。趨勢(shì)分析感染人數(shù)是持續(xù)上升、平臺(tái)期還是下降這對(duì)應(yīng)疫情的發(fā)展階段。異常點(diǎn)檢測(cè)某天的病毒濃度突然飆升或驟降是真實(shí)的疫情暴發(fā)/管控起效還是采樣誤差、檢測(cè)異常、或大型活動(dòng)如音樂(lè)節(jié)、體育賽事導(dǎo)致的人口臨時(shí)聚集造成的我們需要利用統(tǒng)計(jì)方法如3σ原則、移動(dòng)平均法、時(shí)間序列分解識(shí)別異常點(diǎn)并結(jié)合收集到的“上下文信息”如是否舉辦大型活動(dòng)、是否發(fā)布新的防控政策進(jìn)行歸因分析。這一步是將數(shù)據(jù)轉(zhuǎn)化為情報(bào)的關(guān)鍵?？臻g異質(zhì)性分析如果有多點(diǎn)數(shù)據(jù)可以分析不同區(qū)域如商業(yè)區(qū)、居民區(qū)、大學(xué)城的病毒濃度差異。這能幫助定位高風(fēng)險(xiǎn)區(qū)域?qū)崿F(xiàn)精準(zhǔn)防控。2.3 第三層措施評(píng)估——如何量化防控策略的效果這是模型的最終應(yīng)用層。題目可能會(huì)假設(shè)幾種情景例如情景A保持現(xiàn)有措施。情景B實(shí)施全員核酸篩查。情景C對(duì)特定區(qū)域進(jìn)行封控。情景D提倡居家辦公、關(guān)閉娛樂(lè)場(chǎng)所。我們需要利用建立好的模型來(lái)預(yù)測(cè)不同措施下污水病毒濃度及反演的感染人數(shù)隨時(shí)間的變化。這通常需要引入機(jī)制模型如經(jīng)典的傳染病模型SIR, SEIR與我們的反演模型進(jìn)行耦合。例如封控措施降低了人群接觸率這會(huì)在SEIR模型中體現(xiàn)為傳播系數(shù)β的下降進(jìn)而改變每日新增感染人數(shù)I(t)最終映射到模型估算的污水病毒濃度上。通過(guò)比較不同情景下的預(yù)測(cè)曲線如感染峰值、達(dá)峰時(shí)間、總感染人數(shù)可以定量評(píng)估各種措施的成本效益。3. 模型構(gòu)建與求解一個(gè)融合了統(tǒng)計(jì)與機(jī)制的混合框架基于以上拆解我們團(tuán)隊(duì)沒(méi)有采用單一的模型而是設(shè)計(jì)了一個(gè)“分階段、混合式”的建模框架。這個(gè)框架的核心思想是先利用數(shù)據(jù)驅(qū)動(dòng)方法校準(zhǔn)基礎(chǔ)模型再引入機(jī)制模型進(jìn)行推演。3.1 階段一基于貝葉斯推斷的校正因子動(dòng)態(tài)校準(zhǔn)模型面對(duì)反演公式中巨大的參數(shù)不確定性我們選擇了貝葉斯統(tǒng)計(jì)模型作為核心工具。它的優(yōu)勢(shì)在于能將我們對(duì)參數(shù)的不確定性先驗(yàn)知識(shí)和觀測(cè)數(shù)據(jù)的不確定性似然函數(shù)結(jié)合起來(lái)得到參數(shù)完整的概率分布后驗(yàn)分布而不僅僅是一個(gè)點(diǎn)估計(jì)。模型設(shè)定觀測(cè)數(shù)據(jù)時(shí)間序列的污水病毒濃度 (C_t)以及可能不可靠但可作為參考的同期臨床新增報(bào)告病例數(shù) (I_t^{reported})。待估參數(shù)綜合校正因子 (K_t)。我們?cè)试S它隨時(shí)間緩慢變化以捕捉病毒排毒規(guī)律、檢測(cè)效率等隨時(shí)間的變化。我們假設(shè) (\log(K_t)) 服從一個(gè)隨機(jī)游走過(guò)程。狀態(tài)空間模型狀態(tài)方程描述參數(shù)如何演化(\log(K_t) \log(K_{t-1}) \epsilon_t, \quad \epsilon_t \sim N(0, \sigma^2))。觀測(cè)方程描述數(shù)據(jù)如何生成對(duì)于污水?dāng)?shù)據(jù)(C_t (N_t^{infected} / K_t) \nu_t)。這里 (N_t^{infected}) 需要從一個(gè)簡(jiǎn)單的傳播模型如SIR初步估算或者初期可假設(shè)與 (I_t^{reported}) 成比例。(\nu_t) 是觀測(cè)誤差。對(duì)于臨床數(shù)據(jù)如果可用(I_t^{reported} \sim \text{Binomial}(N_t^{infected}, \rho))。其中 (\rho) 是報(bào)告率也是一個(gè)待估參數(shù)。求解與實(shí)現(xiàn) 我們使用馬爾可夫鏈蒙特卡洛MCMC方法通過(guò)Python的PyMC3或Stan庫(kù)來(lái)實(shí)現(xiàn)貝葉斯推斷。代碼的核心是定義上述概率模型然后讓采樣器去探索參數(shù)的后驗(yàn)分布。# 偽代碼示例展示PyMC3模型結(jié)構(gòu)思路 import pymc3 as pm with pm.Model() as wastewater_model: # 先驗(yàn)分布對(duì)log(K)的初始值和波動(dòng)性進(jìn)行假設(shè) sigma_K pm.HalfNormal(sigma_K, sigma0.5) # K隨時(shí)間變化的波動(dòng)大小 K_log pm.GaussianRandomWalk(K_log, sigmasigma_K, shapelen(days)) K pm.Deterministic(K, pm.math.exp(K_log)) # 轉(zhuǎn)換為正數(shù) # 報(bào)告率的先驗(yàn) reporting_rate pm.Beta(reporting_rate, alpha2, beta8) # 假設(shè)報(bào)告率較低先驗(yàn)集中在0.2左右 # 連接感染人數(shù)與污水濃度 (簡(jiǎn)化版假設(shè)感染人數(shù)已知為I_estimated) # I_estimated 可以從一個(gè)獨(dú)立的SEIR模型生成作為輸入 expected_concentration I_estimated / K # 似然函數(shù)假設(shè)觀測(cè)濃度服從正態(tài)分布方差未知 conc_sigma pm.HalfNormal(conc_sigma, sigma10) observed_conc pm.Normal(observed_conc, muexpected_concentration, sigmaconc_sigma, observedreal_wastewater_data) # 連接感染人數(shù)與臨床報(bào)告病例 (如果數(shù)據(jù)可用) reported_cases pm.Binomial(reported_cases, npm.math.ceil(I_estimated), # 感染人數(shù)取整 preporting_rate, observedclinical_reported_data) # 運(yùn)行MCMC采樣 trace pm.sample(2000, tune1000, cores2, return_inferencedataFalse)運(yùn)行這個(gè)模型后我們得到的不是單一的K值而是K隨時(shí)間變化的一系列可能軌跡后驗(yàn)分布樣本。我們可以用其中位數(shù)軌跡作為最佳估計(jì)同時(shí)其分位數(shù)如95%置信區(qū)間清晰地展示了估算的不確定性。這比給出一個(gè)孤零零的數(shù)字要有力得多。3.2 階段二耦合SEIR模型進(jìn)行預(yù)測(cè)與措施模擬在階段一我們可能用一個(gè)簡(jiǎn)單的趨勢(shì)模型來(lái)生成 (I_estimated)。為了進(jìn)行措施評(píng)估我們需要一個(gè)能反映傳染病動(dòng)力學(xué)的機(jī)制模型。我們引入一個(gè)經(jīng)典的SEIR模型[ \begin{aligned} \frac{dS}{dt} -\frac{\beta(t) I S}{N} \ \frac{dE}{dt} \frac{\beta(t) I S}{N} - \sigma E \ \frac{dI}{dt} \sigma E - \gamma I \ \frac{dR}{dt} \gamma I \end{aligned} ]其中(S, E, I, R) 分別代表易感者、潛伏者、感染者、康復(fù)者(N)為總?cè)丝凇?\sigma)是潛伏期倒數(shù)(\gamma)是感染期倒數(shù)。關(guān)鍵參數(shù)是傳播率 (\beta(t))它是一個(gè)可以受防控措施影響的變量。措施如何量化我們將不同的防控措施映射為對(duì) (\beta(t)) 的調(diào)制封控/靜默大幅降低接觸使 (\beta) 降至一個(gè)很低的基礎(chǔ)水平如原來(lái)的20%。社交距離/口罩令中等程度降低傳播效率使 (\beta) 減少一定比例如原來(lái)的40%-60%。疫苗接種通過(guò)降低易感者比例S和可能降低感染者傳染力來(lái)間接影響有效傳播?？梢栽谀Ｐ椭屑尤胍呙绺采w率變量或直接調(diào)整有效接觸率。模型耦合與預(yù)測(cè)流程參數(shù)校準(zhǔn)利用歷史數(shù)據(jù)包括階段一校準(zhǔn)后的感染人數(shù)估算值使用最小二乘法或MCMC來(lái)擬合SEIR模型的參數(shù)主要是初始的 (\beta_0)。情景模擬設(shè)定一個(gè)時(shí)間點(diǎn) (t_{intervention})從此處開(kāi)始將 (\beta(t)) 的值根據(jù)預(yù)設(shè)措施進(jìn)行調(diào)整。運(yùn)行模型數(shù)值求解SEIR微分方程組得到未來(lái)一段時(shí)間內(nèi)每日的感染人數(shù) (I(t))。反演回污水濃度利用階段一得到的最優(yōu)校正因子 (K)或其后驗(yàn)均值通過(guò)公式 (C_{predicted}(t) I(t) / K)預(yù)測(cè)未來(lái)污水病毒濃度的變化趨勢(shì)。效果對(duì)比繪制不同情景下的 (I(t)) 和 (C_{predicted}(t)) 曲線比較峰值高度、達(dá)峰時(shí)間、曲線下面積近似總感染人數(shù)等指標(biāo)。4. 求解全過(guò)程的關(guān)鍵細(xì)節(jié)與避坑指南紙上談兵總是容易但將上述框架實(shí)現(xiàn)并得出可靠結(jié)果過(guò)程中布滿了“坑”。以下是我們?cè)诮忸}和后續(xù)復(fù)盤(pán)時(shí)總結(jié)的幾個(gè)關(guān)鍵細(xì)節(jié)和避坑點(diǎn)。4.1 數(shù)據(jù)預(yù)處理異常值處理與平滑原始的污水監(jiān)測(cè)數(shù)據(jù)噪聲極大。某天濃度奇高可能是因?yàn)椴蓸訒r(shí)正好有一段高濃度污水團(tuán)流過(guò)某天濃度奇低可能是檢測(cè)失敗。直接使用原始數(shù)據(jù)擬合模型會(huì)導(dǎo)致災(zāi)難性后果。我們的做法滑動(dòng)窗口中位數(shù)濾波對(duì)于時(shí)間序列數(shù)據(jù)我們首先采用滑動(dòng)窗口中位數(shù)例如窗口大小為7天進(jìn)行平滑。中位數(shù)比均值對(duì)異常值更不敏感。這能有效濾除短暫的尖峰或低谷保留趨勢(shì)?；谀Ｐ偷漠惓z測(cè)在建立了初步的模型如簡(jiǎn)單的趨勢(shì)模型貝葉斯反演后計(jì)算每個(gè)數(shù)據(jù)點(diǎn)的標(biāo)準(zhǔn)化殘差。將殘差絕對(duì)值大于3倍標(biāo)準(zhǔn)差的數(shù)據(jù)點(diǎn)標(biāo)記為“待審查異常點(diǎn)”。上下文歸因查閱這些異常點(diǎn)對(duì)應(yīng)日期的“上下文信息表”題目若提供或根據(jù)常識(shí)推斷。如果發(fā)現(xiàn)某天有大型集會(huì)則該高點(diǎn)可能是真實(shí)的如果無(wú)特殊事件則很可能是技術(shù)誤差可以考慮用前后數(shù)據(jù)的插值進(jìn)行替換或在模型中賦予該點(diǎn)一個(gè)更大的觀測(cè)誤差方差降低其對(duì)參數(shù)估計(jì)的影響。4.2 校正因子K的不確定性傳遞這是整個(gè)建模中最精妙也最易出錯(cuò)的部分。K不是一個(gè)固定常數(shù)它的不確定性會(huì)直接“傳遞”到感染人數(shù)估算 (N C \times K) 中。在貝葉斯框架下我們得到了K的后驗(yàn)分布那么N的后驗(yàn)分布自然就是 (C \times K) 的分布。在報(bào)告結(jié)果時(shí)必須報(bào)告置信區(qū)間例如“估算今日感染人數(shù)為 1250人95% CI: 800 - 1800”。只報(bào)告點(diǎn)估計(jì)1250人是嚴(yán)重不完整的會(huì)誤導(dǎo)決策者。在將K用于SEIR模型預(yù)測(cè)時(shí)更穩(wěn)妥的做法是進(jìn)行不確定性傳播分析。例如從K的后驗(yàn)分布中抽取1000個(gè)樣本對(duì)每個(gè)樣本K_i運(yùn)行一次SEIR模型得到一條預(yù)測(cè)曲線。最終你會(huì)得到1000條預(yù)測(cè)曲線形成一個(gè)“預(yù)測(cè)帶”。這個(gè)帶子的寬度直觀地展示了由于基礎(chǔ)參數(shù)不確定性導(dǎo)致的預(yù)測(cè)不確定性。這比只用K的均值跑一次模型要科學(xué)得多。4.3 SEIR模型參數(shù)化與過(guò)擬合陷阱SEIR模型看似只有幾個(gè)參數(shù)(\beta, \sigma, \gamma)但在擬合短期、嘈雜的數(shù)據(jù)時(shí)極易過(guò)擬合。特別是初始條件E(0), I(0), R(0)的設(shè)定影響巨大。我們的經(jīng)驗(yàn)固定部分參數(shù)潛伏期 (\sigma) 和感染期 (\gamma) 可以根據(jù)醫(yī)學(xué)研究設(shè)定為固定值如潛伏期5天則 (\sigma 1/5)感染期7天則 (\gamma 1/7)。這減少了待估參數(shù)。謹(jǐn)慎設(shè)定初始感染者I(0)不要簡(jiǎn)單設(shè)為臨床報(bào)告數(shù)?？梢岳梅囱菽Ｐ凸浪愕某跏几腥救藬?shù)作為I(0)的強(qiáng)先驗(yàn)信息。在貝葉斯框架下可以將I(0)也作為一個(gè)待估參數(shù)并為其設(shè)定一個(gè)以反演估算值為中心、有一定寬度的先驗(yàn)分布。使用正則化或先驗(yàn)信息對(duì)傳播率 (\beta) 施加合理的先驗(yàn)分布如對(duì)數(shù)正態(tài)分布防止其擬合出一些物理意義上不合理的極端值。驗(yàn)證如果數(shù)據(jù)量允許采用時(shí)間序列交叉驗(yàn)證。用前80%的數(shù)據(jù)擬合模型預(yù)測(cè)后20%看預(yù)測(cè)效果。避免使用全部數(shù)據(jù)擬合后直接宣稱模型完美。4.4 措施效果評(píng)估的“反事實(shí)”框架評(píng)估措施效果時(shí)最大的挑戰(zhàn)是我們無(wú)法同時(shí)觀測(cè)到“實(shí)施措施”和“不實(shí)施措施”兩個(gè)平行世界的結(jié)果。我們只能觀測(cè)到實(shí)施措施后的實(shí)際數(shù)據(jù)。因此我們的評(píng)估本質(zhì)上是構(gòu)建一個(gè)“反事實(shí)”場(chǎng)景如果沒(méi)有這項(xiàng)措施情況會(huì)怎樣我們的模型校準(zhǔn)后的SEIR就是用來(lái)生成這個(gè)反事實(shí)場(chǎng)景的工具。具體步驟用措施實(shí)施前的所有數(shù)據(jù)將模型參數(shù)校準(zhǔn)到最佳。假設(shè)措施未實(shí)施即 (\beta) 參數(shù)保持措施前的水平不變運(yùn)行模型預(yù)測(cè)措施實(shí)施后的疫情發(fā)展。這是“反事實(shí)”曲線。將反事實(shí)預(yù)測(cè)曲線與實(shí)際觀測(cè)到的數(shù)據(jù)或經(jīng)過(guò)處理后的污水反演感染人數(shù)進(jìn)行對(duì)比。效果量化計(jì)算兩條曲線之間的差異。例如可以計(jì)算在措施實(shí)施期T內(nèi)反事實(shí)曲線預(yù)測(cè)的總感染人數(shù)與實(shí)際估算的總感染人數(shù)之差這個(gè)差值就是模型估計(jì)的“措施避免的感染人數(shù)”。也可以比較峰值降低的百分比。必須強(qiáng)調(diào)這個(gè)評(píng)估的可靠性完全依賴于模型本身的質(zhì)量和我們對(duì)“若無(wú)措施β保持不變”這一假設(shè)的信心。任何模型評(píng)估都必須附帶對(duì)假設(shè)和不確定性的充分討論。5. 論文寫(xiě)作與結(jié)果呈現(xiàn)如何將復(fù)雜計(jì)算轉(zhuǎn)化為清晰故事數(shù)學(xué)建模競(jìng)賽的論文本質(zhì)上是向評(píng)委講述一個(gè)用數(shù)學(xué)工具解決實(shí)際問(wèn)題的完整故事。寫(xiě)作和圖表呈現(xiàn)與模型本身同等重要。5.1 論文結(jié)構(gòu)邏輯我們的論文大致遵循了以下結(jié)構(gòu)但這并非模板而是邏輯的自然流動(dòng)問(wèn)題重述與分析用自己的語(yǔ)言精煉地復(fù)述問(wèn)題并立即展示我們拆解出的“反演-分析-評(píng)估”三層框架圖。讓評(píng)委一眼看懂你的解題思路。模型準(zhǔn)備說(shuō)明數(shù)據(jù)預(yù)處理步驟平滑、異常值處理、列出所有符號(hào)假設(shè)。介紹SEIR模型和貝葉斯反演模型的基本原理并重點(diǎn)闡述為何選擇它們以及它們?cè)诒締?wèn)題中如何連接。模型求解首先展示貝葉斯反演的結(jié)果給出校正因子K的時(shí)間序列圖帶置信區(qū)間以及由此反演的感染人數(shù)估算圖與臨床報(bào)告數(shù)對(duì)比。這張圖是第一個(gè)亮點(diǎn)直觀展示了模型校準(zhǔn)的合理性和不確定性。然后展示SEIR模型的擬合效果用校準(zhǔn)后的感染人數(shù)數(shù)據(jù)擬合SEIR給出擬合曲線圖。最后進(jìn)行措施模擬用不同的β值模擬不同措施繪制未來(lái)感染人數(shù)和污水濃度的預(yù)測(cè)對(duì)比圖。用陰影表示不確定性帶。結(jié)果分析敏感性分析展示關(guān)鍵參數(shù)如先驗(yàn)分布的設(shè)定、報(bào)告率ρ的變化如何影響最終感染人數(shù)估算。這體現(xiàn)了模型的穩(wěn)健性思考。異常點(diǎn)分析列出檢測(cè)到的異常點(diǎn)并結(jié)合假設(shè)的“上下文信息”給出合理解釋。措施評(píng)估表用表格清晰對(duì)比不同措施下的關(guān)鍵指標(biāo)預(yù)測(cè)值峰值人數(shù)、達(dá)峰時(shí)間、避免感染人數(shù)等。模型評(píng)價(jià)與推廣客觀討論模型的優(yōu)點(diǎn)如能量化不確定性、整合多源數(shù)據(jù)、缺點(diǎn)如對(duì)先驗(yàn)信息依賴、假設(shè)簡(jiǎn)化以及未來(lái)改進(jìn)方向如引入空間網(wǎng)絡(luò)模型、考慮病毒變異等。5.2 圖表可視化技巧一圖勝千言在建模論文中尤其如此。圖1數(shù)據(jù)概覽將污水濃度時(shí)間序列和臨床報(bào)告病例數(shù)畫(huà)在同一張圖上雙Y軸第一時(shí)間揭示數(shù)據(jù)的特征和關(guān)聯(lián)。圖2模型原理示意圖繪制一個(gè)流程圖展示“原始數(shù)據(jù) - 預(yù)處理 - 貝葉斯反演 - SEIR模型 - 措施模擬 - 輸出決策支持”的完整邏輯鏈。圖3貝葉斯反演結(jié)果這是核心。用一條深色線表示K或感染人數(shù)的中位數(shù)估計(jì)用淺色帶狀區(qū)域表示其95%置信區(qū)間。務(wù)必在圖上標(biāo)注出識(shí)別出的異常點(diǎn)。圖4SEIR模型擬合與預(yù)測(cè)用散點(diǎn)表示歷史數(shù)據(jù)或反演估算值用實(shí)線表示模型擬合曲線。預(yù)測(cè)部分用虛線延伸并用不同的顏色和線型來(lái)區(qū)分不同措施情景下的預(yù)測(cè)曲線預(yù)測(cè)不確定性用半透明色帶表示。表1措施效果對(duì)比設(shè)計(jì)一個(gè)簡(jiǎn)潔的表格列是不同措施情景行是關(guān)鍵指標(biāo)如傳播率β降低比例、預(yù)測(cè)感染峰值、達(dá)峰時(shí)間、總感染人數(shù)估算、避免感染人數(shù)等。5.3 代碼與可復(fù)現(xiàn)性雖然論文正文不展示大量代碼但在附錄或提交的附件中清晰、注釋良好的代碼是加分項(xiàng)。我們當(dāng)時(shí)使用了Python主要庫(kù)包括pandas數(shù)據(jù)處理、numpy數(shù)值計(jì)算、matplotlib/seaborn繪圖、pymc3貝葉斯建模、scipy微分方程求解。代碼文件按功能模塊組織并有一個(gè)主腳本清晰地按順序調(diào)用各個(gè)模塊確保評(píng)委或任何人能夠一鍵復(fù)現(xiàn)我們的主要結(jié)果。回顧整個(gè)解題過(guò)程從最初面對(duì)“污水?dāng)?shù)據(jù)”和“校正因子”這些陌生概念的茫然到最終構(gòu)建出一個(gè)能自圓其說(shuō)、量化不確定性的混合模型最大的收獲不是學(xué)會(huì)了某個(gè)特定算法而是鍛煉了一種“結(jié)構(gòu)化定義問(wèn)題”和“用概率思維擁抱不確定性”的能力。在實(shí)際的公共衛(wèi)生決策中數(shù)據(jù)從來(lái)都是不完美的模型也永遠(yuǎn)是現(xiàn)實(shí)的簡(jiǎn)化。一個(gè)優(yōu)秀的模型其價(jià)值不在于做出精準(zhǔn)無(wú)比的預(yù)測(cè)而在于能清晰揭示關(guān)鍵驅(qū)動(dòng)因素、量化不同選擇的風(fēng)險(xiǎn)與收益并將結(jié)論的不確定性坦誠(chéng)地呈現(xiàn)給決策者。這次數(shù)學(xué)建模的經(jīng)歷正是對(duì)這一理念的一次深刻演練。