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