研究筆記

只見總數,能否框住看不見的疫情?

延遲的地區總數可以支持看似精確的估計,卻未必緊緊限制隱藏狀態。透過六區 SIR 合成基準,比較 ensemble 區間與經驗證集合界限,保留界限闊得幾乎失去作用的結果。

2026年9月12日 · 約 25 分鐘閱讀

一個監測畫面顯示幾個鄰近地區的感染人數總和。數字比現在慢了幾天,量測本身也不完全準確。在同一個總數背後,某一區可能已經回落,另一區卻仍在上升。估計算法替六個地區畫出平滑曲線,每條曲線旁邊還有一條窄窄的陰影帶。畫面看起來很具體,也很有把握。

可是,究竟是哪一件事變得具體了?資料真的已經排除了所有差異很大的地區軌跡,還是算法在自己的機率模型之下,選出了一小批特別合理的軌跡?前者談的是還有哪些情況不可能被排除,後者談的是哪些情況比較值得相信。它們可以互相補充,卻不是同一種承諾。

本文用一個假想的六區疫情,研究這個差別。這裏沒有病人紀錄,沒有從真實疾病估計出來的參數,也沒有對當前疫情作出判斷。每一筆觀測都來自已知的合成模型。正因為我們知道被隱藏的真值,才可以逐項檢查:兩種看似相似的不確定性區間,究竟有沒有包含真值,又為了包含它而付出了甚麼代價。

第一種方法是 ensemble Kalman filter,簡稱 EnKF,中文常稱集合卡爾曼濾波。它維持一批模擬狀態,利用觀測量與隱藏變量之間的共同變動來更新估計。第二種是 set-membership estimation,本文稱為集合成員估計:只要某個狀態沒有違反已知初值、模型及誤差界限,就不應把它排除。當精確可行集合太複雜,便用可計算的外包絡保留它。

關鍵在於「外包絡」。一個界限可以成功保留真值,同時也保留太多其他可能性。證明包含真值,不等於證明答案有用;相反,一條很窄的 ensemble 區間可以有實際參考價值,卻不因此成為確定性界限。研究項目簡介提供摘要;下文會沿着模型、八幅研究圖和失敗案例,逐步拆開這些差別。

這個比較發現了甚麼

兩組、延遲五天的基準,最能顯示這個張力。50 次重複中,集合外包絡的平均寬度跨越地區人口的 48.40%,EnKF 則為 0.951%。集合方法保留了每一條完整感染軌跡,EnKF 沒有一條能在全部日期、全部地區都不漏掉真值,但其平均分量包含率仍有 94.05%。這不是互相矛盾的數字,而是不同的可靠性問題。

第 20 日發出預測後停止收取新資料,基準集合峰值大小範圍會變成 0–6,000,候選日期跨越 21–120。包含性仍在,精度卻幾乎消失。因此,本文不會把某個方法封為普遍勝者,而是分開討論分布資訊的用途、條件式保護,以及算法本身的保守性。它們需要不同的檢查,也可能適合不同問題。

總數是一項量測,不是一張地區分布圖

這個實驗觀察的是「當時仍然感染的人數」,不是「一天之內新增的感染」。這個分別決定量測方程應該怎樣寫。感染存量描述某一時刻的狀態;每日新增病例則計算一段時間內有多少人進入感染狀態。即使量測完全沒有誤差,康復、人口流動和報告時間也會令兩個數量不同,不能因為兩者都以人數表示,就互相代替。

六個地區的初始人口依次為 1,000、900、1,100、950、1,050 和 1,000。初始感染人數的中心依次是 8、5、6、4、7 和 5,每區可獨立地比中心少兩人或多兩人。初始康復人數全部為零,其餘人口屬於易感者。往後模擬出現小數人數,是連續分艙微分方程的近似,不代表現實存在「半個感染者」。

我們比較三種報告方式。六組觀測逐區報告感染存量;兩組觀測分別報告第一至第三區,以及第四至第六區的總和;一組觀測只報告六區總數。每一組的量測都是所屬地區感染存量的總和,再加上一個有已知絕對上限的誤差。

六區合成感染真值曲線,以及合成兩組和一組之後按收到日期顯示的延遲觀測,說明總數保留和丟失的地區資訊。
圖 1.預先選定的 replicate 0 中,六區 SIR 實驗究竟觀察到甚麼。上圖分開六個感染真值,下兩圖把它們合成兩組或一組。兩組面板中,藍色圓點代表地區 1–3,橙色方塊代表地區 4–6。曲線是當時真值,空心標記是在收到當日顯示、延遲五天且每組誤差不超過五人的報告。為方便閱讀,標記只顯示每第五份報告;算法使用全部每日報告。基準過程上限為 0.5,每日流動率為 0.002。

假設第一組有 19 名感染者。分布為 (8, 5, 6) 或 (10, 3, 6),總數都一樣。這兩個向量都落在實驗已聲明的初始不確定範圍之內。第二個分布不是濾波器的輸出,也不是為了展示成功而挑出來的估計,而是一個用來說明量測映射的可行初始狀態。其易感人數作相反調整,因而保持每區的初始人口不變。

單看這個總和,資料沒有保留「第一區多兩人、第二區少兩人」的方向。只要組內此消彼長,總數就察覺不到差異。如果我們永遠只能知道這一個數字,這種歧義便無法消除。不過,本文供應的是一段時間序列,而且各區感染率不同;兩種初始分布往後的總數,可能因動力學不同而逐漸分開。

因此,單次聚合量測有歧義,不等於整個動力系統不可觀測。 某一刻沒有保留的資訊,可能在日後的變化中留下痕跡。我們沒有證明六區模型存在全域的不可識別性,也沒有聲稱只要聚合便一定無法重建地區狀態。基準測試回答的是較具體的問題:兩個已固定的算法,如何利用一段有噪聲、又有延遲的聚合觀測。

兩個符合初值範圍的感染分布,在第一及第二區不同,但兩組報告的總數完全重合,展示單次聚合量測的歧義。
圖 2.精確的單次量測反例,不是 Monte Carlo 結果。狀態 A 的感染人數為 (8,5,6,4,7,5),狀態 B 把第二區兩名感染者轉到第一區,並以相反的易感調整保持各區人口,兩者均符合初始界限。下圖兩組總數完全相同,空心方形與實心圓形刻意重合;上圖的小幅水平位移只為分開類別標記。這並不推出未來整段觀測序列也相同。

比較還有一項容易忽略的細節。誤差上限是附在每一組報告上,不是附在每一個人或每一個地區上。一組總數有一個有界誤差;把六份分區報告相加,則會把六個誤差合起來。因此,三種方式是不同的資訊系統,並非已經控制感測器成本的採購實驗。改變分組既改變空間解析度,也改變總數的誤差結構,不能把結果直接翻譯成應該購買多少感測器的建議。

人會流動,感染項的分母也會改變

每區包含易感、感染及康復三個狀態,分別寫作 SiS_iIiI_iRiR_i。當時人口是它們的總和 Ni=Si+Ii+RiN_i=S_i+I_i+R_i。六區排成一個環,每人以相同速率流向左右鄰區。第一區的鄰居是第二區和第六區。三個分艙都使用同一個流動規則,並非只有感染者會移動。

在每日擾動之間,模型為

Mi(Z)=m(Zi1+Zi+12Zi),S˙i=βiSiIiNi+Mi(S),I˙i=βiSiIiNiγIi+Mi(I),R˙i=γIi+Mi(R).\begin{aligned} \mathcal M_i(Z)&=m(Z_{i-1}+Z_{i+1}-2Z_i),\\ \dot S_i&=-\beta_i\frac{S_iI_i}{N_i}+\mathcal M_i(S),\\ \dot I_i&=\beta_i\frac{S_iI_i}{N_i}-\gamma I_i+\mathcal M_i(I),\\ \dot R_i&=\gamma I_i+\mathcal M_i(R). \end{aligned}

感染率依次為每天 (0.36,0.34,0.38,0.33,0.35,0.37)(0.36,0.34,0.38,0.33,0.35,0.37),康復率為每天 0.140.14。基準設定中,流向每一個鄰區的人均速率為每天 0.0020.002。這些都是合成模型的固定輸入,不是任何疾病的實測估計。濾波過程不會再辨識未知參數,亦沒有以更換參數來配合不理想的結果。

把同一區的三條方程相加,感染和康復項會抵消,只剩下人口流動:

N˙i=m(Ni1+Ni+12Ni),i=16Ni(t)=6000.\begin{aligned} \dot N_i&=m(N_{i-1}+N_{i+1}-2N_i),\\ \sum_{i=1}^{6}N_i(t)&=6000. \end{aligned}

守恆的是六區總人口,而不是每一區的人口。若一直使用初始人口作感染項的分母,就等於悄悄換了模型。在這個對稱環上,流動具有平均化作用:人口最多的區不能再向上增加,人口最少的區不能再向下減少。因此,各區人口一直介乎 900 和 1,100 之間。這項物理約束也令分母遠離零,為後面的包含論證提供必要條件。

每完成一天的連續演化,我們會在同區易感者與感染者之間施加一次有界轉移。只有當提出的轉移會令分艙變成負數,才把轉移截到可行範圍內。基準設定中,每區每天最多轉移半人,可以增加或減少感染。它是刻意簡化的模型差異項,不應解讀成某一種特定流行病學機制的寫實描述。

ui,n=clip(u~i,n,Ii(n),Si(n)),Si(n)=Si(n)ui,n,Ii(n)=Ii(n)+ui,n,Ri(n)=Ri(n).\begin{aligned} u_{i,n}&=\operatorname{clip}\bigl(\widetilde u_{i,n},-I_i(n^-),S_i(n^-)\bigr),\\ S_i(n)&=S_i(n^-)-u_{i,n},\\ I_i(n)&=I_i(n^-)+u_{i,n},\\ R_i(n)&=R_i(n^-). \end{aligned}

其中 nn^- 表示每日轉移發生前。提出的轉移滿足 u~i,nq|\widetilde u_{i,n}|\le q,基準為 q=0.5q=0.5。易感者減少多少,感染者就增加多少,反方向亦然,故每區人口在這一步不變。然後才量測轉移後的感染存量,再經過報告延遲送到估計器。次序是連續演化、有界轉移、量測,以及延遲送達。

把這個次序寫清楚,可以避免幾種看似細小、其實改變了問題的做法。如果給三個分艙各自加入獨立誤差,人口便不再守恆;把每日有界轉移當成連續白噪聲,會換掉不確定性模型;如果每個數值子步都允許一次完整擾動,求解器步長愈小,整天允許的擾動總量便愈大。本文沒有採用這些替代做法。

靠近非負邊界時,截斷也會改變實際擾動的分布。感染人數接近零,負向轉移不能移走多於現有的感染者,正向轉移卻仍可把易感者移入感染分艙。所以,提出的誤差雖然對稱,實際執行後未必仍然是零均值。真值與 ensemble 使用同一規則,集合方法則包含所有允許的截斷轉移。若這個機制支持了模擬後段的小量感染,也不能把它當成已觀察到再感染或外來輸入病例的證據。

一份報告有兩個日期

每份報告都有「量測所屬日」和「可以使用的日子」。延遲五天時,第 20 日收到的報告量測的是第 15 日。誤差上限約束第 15 日的總數,不能因為報告今天才到,就把它的量測條帶直接套在第 20 日的狀態上。

若分組矩陣為 CC,量測方程為

yn=CI(n)+vn,vg,nvˉ,arrival day=n+d.\begin{aligned} y_n&=C I(n)+v_n,\qquad |v_{g,n}|\le\bar v,\\ \text{arrival day}&=n+d. \end{aligned}

矩陣每一列在所屬地區填一,其餘位置填零。延遲 dd 已知,且只接受整日;實驗比較 0、5 和 10 天。非整日輸入會被拒絕,不會靜默向下取整到某一天。研究開始前沒有虛構的歷史量測,第一份報告到達之前,估計器只能靠初始資訊和模型推進。

兩個算法都把已同化的狀態保存在最新的量測所屬日。下一份報告到達時,先把這個歷史狀態推進一天,在正確日期加入量測,再跨越已知延遲推進到現在。最後這一段只使用模型和允許的擾動,不會讀取明天才收到的報告。兩個方法共享同一批數據,也共享同一條資訊時間線。

這不只是資料整理習慣。如果疫情正在上升,把較舊、較低的總數當成當前量測,可能把今天的估計拉低;疫情正在下降時,同一種錯誤可能把估計拉高。這樣得到的窄帶不代表資料特別精確,而是時間對齊錯誤被畫成精確結果。比較必須先排除這種假象,才有資格討論算法本身的表現。

延遲也改變運算量。每次在歷史日期更新之後,十日延遲需要重新向前推進十天,才能回答今天的問題。歷史後驗狀態與當前預測並不是同一件東西。因此,報告速度時必須說明是否包含這段額外推進。本文的濾波計時包括它,但不包括另外進行的長期峰值預測。

EnKF 從成員之間學到甚麼

EnKF 先維持許多可能狀態,主實驗使用 500 個 members。每個成員從與真值相同的初始感染界限內獨立抽樣,之後接受自己的一系列有界守恆擾動。成員之間逐漸形成的分布,反映了初始假設與過程模型如何把不確定性帶到未來。

收到報告時,每個成員都先預測對應的地區總數。假如總數較大的成員通常也有較高的第二區感染人數,這種共變動便讓濾波器可以透過一個總數去更新第二區。這是有用的結構,但其來源部分是 ensemble 的分布與模型演化,不全部來自剛收到的量測。理解這一點,才不會把算法給出的六個估計誤認為六項獨立觀測。

分析更新寫成

K=Pxy(Pyy+Rv)1,xa,(j)=xf,(j)+K(y+ε(j)CIf,(j)),Rv=vˉ23Id.\begin{aligned} K&=P_{xy}(P_{yy}+R_v)^{-1},\\ x^{a,(j)}&=x^{f,(j)}\\ &\quad+K\bigl(y+\varepsilon^{(j)}-C I^{f,(j)}\bigr),\\ R_v&=\frac{\bar v^2}{3}\,\mathrm{Id}. \end{aligned}

PxyP_{xy}PyyP_{yy} 是成員之間的樣本協方差,上標分別表示預測與分析狀態。真實量測誤差及更新時的人造量測擾動,均在 [vˉ,vˉ][-\bar v,\bar v] 上作均勻抽樣,各組相互獨立。這個分布的變異數是 vˉ2/3\bar v^2/3;若直接把誤差上限平方當成變異數,便是在測試另一個濾波器。初始誤差與提出的每日擾動,同樣從已聲明的範圍內均勻抽樣。

線性更新有可能產生負數或違反人口約束,因此每次分析後,都會把每個成員在每區的三個分艙,投影到非負、且總和等於當時區內人口的單純形上。這裏的投影是歐氏距離下最近的可行狀態,人口由獨立的線性人口方程推進。這是預先聲明的算法步驟,不是事後刪去不好看的成員。

本次沒有協方差膨脹,也沒有局部化。這讓比較保持固定而容易檢查,同時也限制結論:我們研究的是這一個 EnKF 設定,而不是宣稱已找到調校得最好的 ensemble 方法。改變成員數可以檢查抽樣敏感度,卻不能代替對所有校準技巧的完整研究。

圖中 EnKF 的上下界,是每個感染分量在成員之間的第 2.5 和第 97.5 百分位,即中央 95% ensemble 區間。在這個非線性、又有物理投影的更新程序裏,它們不是精確貝葉斯可信區間,更不是所有符合誤差界限狀態的確定性包絡。即使是一個理想的中央機率區間,也刻意把部分機率留在兩端之外。

此外,單一分量的區間與整條軌跡的區間承諾也不同。第 15 日第二區落在帶內,不保證第 38 日第四區也落在帶內。要求每一天、每一區都不漏掉真值,遠比平均包含率接近 95% 嚴格。我們會同時報告這兩種指標,而不把每次模擬中的 720 個地區日檢查當作互相獨立的實驗。

集合界限拒絕丟掉甚麼

集合成員估計問的是:如果模型及所有誤差上限都正確,哪些狀態仍然有可能?它不必指定界限之內的機率密度。靠近邊緣的擾動與靠近中心的擾動,同樣不能無故排除。用來生成本次合成資料的均勻分布,並不是集合包含論證的必要條件;論證依賴的是支持範圍。

精確可行集合很快會變得複雜。非線性動力學把它彎曲,反覆量測把它切開,人口流動又把地區互相連接。這裏使用較簡單的表示:替每個分艙保留下端點與上端點,形成箱形外包絡,並帶着非負與人口守恆約束。這是刻意保守的基準,不是宣稱箱形已是最有效率的集合表示。

如果一組感染總數介乎 aabb,其中某一區的感染人數便必須滿足

IjaijUi,IjbijLi.I_j\ge a-\sum_{i\ne j}U_i,\qquad I_j\le b-\sum_{i\ne j}L_i.

這裏的總和只包含同組其他地區。把新上下限與原本區間取交集,不會刪去符合群組量測的狀態。同一原理也用於 Si+Ii+Ri=NiS_i+I_i+R_i=N_i 和總人口 6,000 的約束。反覆收縮可能得到更好的端點,但固定次數的收縮未必找到最緊的箱形外包絡。

這個簡單規則也解釋了總數的限制。知道一組有 19 人,究竟對第一區有多少幫助,取決於另外兩區還可能貢獻多少。其他地區的不確定範圍愈闊,就愈能透過組內重新分配保持同一總數。界限所表達的是還不能排除哪些替代狀態,不是說每種狀態一樣可能。

量測之間的傳播更加需要小心。只積分箱形的下角點和上角點,通常不能保證夾住非線性疫情模型;變量之間的耦合可能有相反方向的效果。抽樣大量軌跡也不足以證明沒有任何未抽中的可行軌跡逃出包絡。這正是「數值上見過很多例子」與「對所有允許情況的包含論證」之間的分界。

算法會先在目前箱形 BB 周圍建立一條管狀區域 YY,再用區間運算計算其中物理可行部分的向量場範圍 F(Y)F(Y)。只有在步長 hh 滿足以下嚴格包含檢查時,才接受這一步:

B+[0,h]F(Y)int(Y),BnextB+hF(Y).\begin{aligned} B+[0,h]F(Y)&\subset\operatorname{int}(Y),\\ B_{\mathrm{next}}&\subseteq B+hF(Y). \end{aligned}

第二個式子的意思,是下一個箱形由向外捨入的端點包絡,再與已知物理約束相交而來;包含性談的是符合約束的真實軌跡,不是被約束刪掉的任意箱內點。第一個式子提供關鍵論證:假設物理軌跡第一次離開管狀區域,在離開之前,它的導數都已包含在所計算的範圍內,累積變化仍會把它放在管狀區域嚴格內部,與「第一次離開」矛盾。

因此,端點包絡保留了這一步內軌跡可能經歷的所有導數積分。它是一階的範圍包絡,不是漏算截斷誤差的普通 Euler 步進。代價就是闊度:把隨時間改變的導數替換成一個整段範圍,會失去時間與變量之間的關係;較高階的經驗證方法可能保留更多這類資訊。

區間加、減、乘、除都在有限 IEEE 754 雙精度運算假設下逐次向外捨入,總和亦使用明確捨入的加法,並檢查分母為正。一般步長為四分之一天;如果管狀包含測試失敗,就細分步長。到達最小允許步長仍失敗,便回報驗證錯誤,不能把未驗證的數值當成界限繼續使用。

每日擾動以容許增量範圍加到易感與感染分艙,再用物理約束收縮。這個外包絡丟失了「兩個增量必定相反」的關係,因而可能變闊,卻不能以這種資訊丟失為理由把它收窄。這些相依關係在每次傳播中反覆丟失,是後面保守界限的重要來源。

完整論證必須逐步跟隨實際運算:初始集合包含允許起點;經驗證傳播保留包含性;有界守恆轉移保留它;在正確歷史日期加入量測條帶保留它;跨越延遲推進到現在亦保留它。如果交集為空,就保留空集合結果,不會偷偷放大量測誤差,或從已知真值重新開始。

怎樣量度有用,而不是把同一個實驗重複計數

主設計交叉比較三種分組、三種延遲及兩種量測誤差上限:6、2 或 1 組;0、5 或 10 天;每份報告最多誤差 1 或 5 人。共 18 個設定,每個設定有 50 次獨立重複。另以兩組、延遲五天、誤差五人為基準,單獨把過程擾動上限改成 0 或 1,或把流動率改成 0 或 0.01,加入四個敏感度設定。合計 1,100 次完整的 120 日模擬重複。

每次重複中,兩個方法使用配對資料。不同亂數流負責初始真值、過程擾動、量測及 ensemble;相容的設定共用底層隨機抽樣,避免比較無謂地被不同起始疫情主導。統計上的獨立單位仍是整次 replicate,不是其中某一天、某一區,亦不是濾波器內某一個成員。

分量包含率先在每次重複內,對六個感染狀態及第 1 至第 120 日的包含指示取平均。整條軌跡同時包含率則只有在全部檢查都成功時才記一。區間寬度先除以當天實際地區人口,然後才平均:

w=1720n=1120i=16Ui(n)Li(n)Ni(n).w=\frac{1}{720}\sum_{n=1}^{120}\sum_{i=1}^{6}\frac{U_i(n)-L_i(n)}{N_i(n)}.

正規化寬度為 0.5,表示平均一個分量的區間跨越了當時區內人口的一半,不是說感染機率為 50%。平均寬度也不能代表每一天:直接觀測剛更新時可能很窄,長時間預測之後卻可能幾乎沒有限制。讀圖時要同時看整段軌跡和彙總指標,不能只選有利的一種。

平均指標的區間使用 2,000 次整個 replicate 的百分位 bootstrap 重抽樣;方法差異採配對重抽樣。這樣保留每次模擬內的時間與地區相依性。如果 50 次全部成功,普通經驗 bootstrap 的上下限會同樣等於一。這是重抽樣計算的性質,不是證明母體成功機率恰好為一。集合方法若有確定性的條件式承諾,其根據是前述包含論證,而非這個退化的統計區間。

運算時間以一批向量化模擬的實際經過時間,除以 replicate 數目及更新次數,當中包含由歷史更新推進到現在的工作。這是本機實作的攤分吞吐量,不是單一即時個案的延遲。整批只共享一個計時結果,所以不會把它複製成 50 個「獨立時間樣本」,再畫出沒有根據的時間信賴區間。

參考真值由高精度自適應積分產生,並與更嚴格容差的結果比較。Ensemble 使用四階 Runge–Kutta 及四分之一天步長,再與一半步長核對。這些收斂檢查處理的是數值近似問題,不會把 ensemble 百分位升格為保證,也不是用真實疾病資料驗證模型。

隱藏狀態的區間實際長甚麼樣

逐區、即時的報告能令簡單集合表示保持相對緊。固定每份量測的誤差上限為五人,只增加延遲,便得到下列寬度。所有數值都是相對於當時地區人口的百分比,再對 50 次重複取平均:

六組報告的延遲集合寬度 (%)EnKF 寬度 (%)
0 天0.7170.268
5 天8.4240.378
10 天26.5270.536

集合方法的急劇增闊,不全是物理疫情突然多了不確定性。每一個沒有量測的步驟,也會令一階箱形表示丟失一些相依關係。在五日延遲下,把六份分區報告改成兩組總數,平均寬度更升至 48.40%。圖 3 顯示預先選定的一條軌跡如何出現這種差異;圖 4–5 則保留較完整的設定與重複比較,不能只靠一幅例圖概括所有情境。

同時,整條軌跡失敗,也不應被誤解成分量估計完全沒有用。上表三個六組設定的 EnKF 分量包含率,依次約為 94.53%、95.09% 和 95.21%;可是,每個設定的 50 條軌跡,都沒有一條能在每一天、每一區全部落在區間內。平均包含率甚至不是一直低於 95%。若問題需要同時保護整條軌跡,就要建立對應的同時不確定性構造,不能替原本的邊際百分位換個名稱便當成做到了。

十次配對重複的 ensemble-size 檢查中,兩組基準使用 100、500 和 1,000 個成員,分量包含率分別為 90.03%、92.67% 和 92.50%。增加成員可以減少其中一種抽樣不確定性,卻不保證包含率單調改善。這十次是敏感度子集,不是拿來取代主要 50 次估計的較小樣本;也不能因為某個成員數在子集較好,就把它事後改成主設定。

第一區的真值、橙色集合外包絡與藍色中央 ensemble 區間,比較逐區即時報告,以及兩組聚合並延遲五天的情況。
圖 3.預先選定 replicate 0 的第一區,量測誤差上限五人、過程上限 0.5、每日流動率 0.002。上圖有六份無延遲分區報告,下圖只有延遲五天的兩組總數。黑線是合成真值,橙色虛線是受約束外集合的界限,藍色實線是 500 個成員的中央 95% 區間。陰影不是平均表現的信賴區間。兩圖使用不同縱軸尺度,避免其中一個情境被壓成扁平線條。
每種方法全部九百個主要重複的包含率與寬度,分開顯示分量包含,以及每一區每一天均包含的較嚴格整條軌跡事件。
圖 4.每種方法全部 900 個主要 replicate 摘要,即 18 個設定各重複 50 次。橫軸以對數尺度顯示平均區間寬度除以當時地區人口;上圖是分量包含率,下圖是第 1–120 日的整條軌跡包含指示。方形代表集合方法,圓形代表 500-member EnKF。沒有為了避免重疊而移動數值,尤其零或一的軌跡分數會刻意重合。地區日與成員都沒有被當成獨立重複。
四個附數字矩陣分開比較兩種方法的平均人口正規化寬度與分量包含率,橫跨三種分組及三種報告延遲。
圖 5.量測誤差上限五人、基準過程上限 0.5 和每日流動率 0.002 時,分組與延遲的比較。每格是 50 次獨立重複的平均;首兩圖顯示正規化寬度,後兩圖顯示分量包含率。黑色數字承擔定量比較,淺色背景只協助閱讀。包含率 1.000 是本次有限實驗的紀錄,不是普遍成功的統計證明。誤差上限一人的設定仍保留在圖 4 及完整比較中。
攤分濾波時間,以及每日過程上限和人口流動的單因素敏感度,將區間寬度與包含率分開並以整次重複估計統計區間。
圖 6.上圖:兩組、量測誤差五人時觀察到的濾波吞吐量,包含跨越報告延遲的推進,合併十個一批的計時,不代表單一個案延遲,亦不虛構時間誤差棒。下兩圖:基準與單獨改變 q=0、1 或 m=0、0.01 的設定,全部延遲五天、各 50 次重複。上方與中間面板的縱軸採對數尺度。誤差棒是 2,000 次 replicate 重抽樣的 95% 百分位 bootstrap 區間。改變 q 同時改變合成擾動與允許的不確定性,不只是調校估計器。

即使關掉過程擾動,也未能挽救這個箱形表示。在同樣的聚合、延遲設定下,把 q 設為零,平均正規化寬度仍有 47.64%,只比基準的 48.40% 略低。初值不確定性、延遲總數及丟失的相依關係,已足以令界限保持寬闊。這個比較同時改變容許的擾動和生成的真值,並不是測試估計器錯誤忽略真實擾動時會怎樣。

寬度代價也不同於時間代價。這次基準運算記錄的攤分吞吐時間,集合方法約為每次重複、每次更新 5.0 ms,EnKF 則約為 15.7 ms。在這個實作中,較窄的輸出反而用了較多運算。這些數字來自當次工作負載下的批次計時,不是獨立個案延遲,也不是不受硬件影響的排名。更精細的集合表示可能同時改變寬度與成本。

未來峰值是另一個問題

第 20 日時,我們固定兩個估計器當時可用的資訊,發出一次預測,此後不再同化新量測。若延遲十天,最新收到的報告只量測到第 10 日。這樣得到的是當時真正能做的預測,不是事後看過峰值後再平滑回去的估計。

目標是第 21 至第 120 日每日網格上的總感染峰值。如果每天總感染人數的下、上包絡為 LnL_nUnU_n,峰值大小與候選峰值日由以下關係得到:

maxnLnmaxniIi(n)maxnUn,D={n:UnmaxsLs}.\begin{gathered} \max_n L_n\le\max_n\sum_i I_i(n)\le\max_n U_n,\\ \mathcal D=\{n:U_n\ge\max_s L_s\}. \end{gathered}

真正最大值所在的日期一定屬於 D\mathcal D,因為該日上界不能低於另一日已經被下界強迫達到的數值。候選日期可以不連續。為了圖像簡潔,會顯示最早至最遲的外包範圍,但包含率使用保存下來的完整候選日期集合檢查。外包範圍可能包含原集合已排除的中間日子,不代表模型預測了一整段平坦高峰。

總感染包絡也保留全域人口守恆,不能超過 6,000;易感及康復人數的下界,可以進一步減少感染上限。只把六區各自的上端點相加,會忘記這六個最大值未必可以同時出現。把已知總和約束用於推導的總感染範圍,與因為圖太闊而人為剪裁,是兩件不同的事。

EnKF 則先讓每個成員產生完整的未來軌跡,各自計算總感染峰值及達峰日,之後才在成員之間取中央百分位。如果先對每天取百分位,再找百分位曲線的最大值,回答的是另一個問題。所有日期都只指每日網格,並不聲稱已限制兩天之間的最大值,或得到連續時間的精確峰值日期。

百分位插值有時會把 ensemble 日期端點放在兩個整數之間。那只是整數日期分布的數值摘要,不是多出了不足一天的時間資訊。同樣地,很多成員可以在同一天達峰,令日期區間寬度為零,卻仍對峰值大小有不同預測。這種表面上的日期精確度,需要用之後的真值檢查,不能只因區間很窄便當成可靠。

兩組基準的 50 次集合預測,全部保留 0–6,000 的峰值大小區間及 21–120 的日期跨度。EnKF 的峰值大小區間在 49 次中包含真值,日期區間在 50 次中全部包含真值,平均寬度約為 53.06 人及 0.44 天。replicate 0 的真正每日網格峰值為 1,424.32 人,出現在第 26 日;其 ensemble 日期區間同樣是第 26 日。

這個有利的結果屬於已固定的合成模型,不是對一般疫情預測的承諾。這裏已知感染率、康復率、流動和延遲,未來擾動也遵守生成資料時的規則。如果現實中其中一項條件改變,不能只搬用這個包含率。相反,集合範圍雖然在這個例子太闊,亦不能由此推斷所有經驗證方法都無法限制峰值。

第20日發出的預測,顯示闊大的總感染外包絡、五十次重複的 ensemble 峰值大小區間誤差,以及差異很大的每日峰值日期範圍。
圖 7.兩組、延遲五天、量測誤差五人及基準動力學。預測在第 20 日發出,只使用當時收到的報告,此後不再同化量測。上圖:replicate 0 在第 21–120 日的真值與集合總感染包絡。中圖:全部 50 次 EnKF 峰值大小區間,減去各自真實峰值;沒有跨過零的線段就是漏掉真值。下圖:replicate 0 的峰值日期範圍,點線標出真正每日網格峰值日。集合日期範圍是保留候選日期的外包區間。

界限可能錯、可能空,也可能只是沒有用

三種失敗值得分開命名。第一種是誤差假設不正確,卻仍得到非空集合。下面的單步例子中,第一區真正有八名感染者;實際誤差 0.5、聲明上限一人,得到 [7.5,9.5]。若實際誤差是 1.5,仍用一人的上限分析,便得到 [8.5,10]。答案看起來有條理,集合也沒有變空,卻已排除了八這個真值。非空本身並不驗證假設。

第二種是互相衝突的報告令交集變空。在已收縮的狀態上,再加入偏移十人的報告,便產生這個結果。實作會回報矛盾,不會修改噪聲容許量。空集合表示採用的約束不能同時成立,但它本身不能指出究竟是哪個感測器、哪個模型項,或哪項初始假設出了問題。直接選一個最順眼的答案,反而會丟失這個診斷訊息。

第三種是運算仍然有效,表示方式卻變得沒有用。補充的解析對照由零感染、零康復及零過程擾動開始,只有易感人口流動,真正感染人數一直恰好為零。不加入新觀測,帶有留白的簡單一階箱形仍會逐步擴大,直到最大的地區感染上界約為 1,063.63。管狀構造容許的微小正範圍被反覆傳播,相依關係又逐步丟失,最終形成非常保守的答案。

這是實作包絡的限制,不是無病狀態本身不可知道的證據。我們其實已經有感染恆等於零的解析解。若把這個特例中的寬界限全部歸因於「資料不足」,便會把算法表示上的損失與真正的資訊限制混在一起。更好的集合表示或特定不變量收縮是否能改善它,是下一個可以單獨測試的問題,不應在看到負面結果後偷偷換方法,再宣稱原方法從未失敗。

真值為八名感染者時,有效與超出假設的量測界限、明確空集合結果,以及完全沒有感染的解析軌跡外逐漸擴大的箱形上界。
圖 8.獨立診斷例子,不計作額外主要 replicate。上圖:第一區真值八人、聲明誤差上限一人,實際誤差 0.5 和 1.5 產生不同保留區間。再加入偏移十人的報告會得到空集合,故該列不畫區間。下圖:q=0、基準流動、120 日沒有新觀測的精確無病對照。黑色零線是真值,虛線是最大的地區箱形上界。這裏沒有附加信賴區間,也不作真實疫情解讀。

文獻如何改變我們應該下的結論

疫情的區間估計不是新問題。Aronna 與 Bliman 已在指定模型及量測假設下建立 SIR–SI 區間觀測器;Gouzé 等人的工作更早研究了不確定生物系統的區間觀測器。因此,不能把「替隱藏疫情建立界限」這個一般概念當成本文的新貢獻。不過,它們的構造與假設,也不等同這裏的混合六區箱形基準。SIR–SI 區間觀測器生物系統區間觀測器

延遲量測同樣已有研究。Meslem 等人討論離散延遲量測下的集合值估計,Efimov 等人分析延遲量測的區間觀測器。這些工作提醒我們必須交代資訊時鐘及包絡成立的假設,卻不會自動替一個已知延遲的非線性實作提供正確性。自己的傳播與收縮仍需要逐步論證。延遲集合值估計延遲相依區間觀測器

表示方法方面,Rego 等人的受約束 zonotope 研究說明:需要保留不變量時,箱形不是唯一選擇。本文箱形很闊,不能證明所有集合方法都一樣沒有資訊。經驗證常微分方程積分也有完整研究脈絡,Nedialkov 等人的工作是其中參考;本文透明的一階包絡,只是這個更大領域裏的一個簡單基準。保留非線性不變量的集合估計經驗證初值問題積分

資料同化方面,Wang 與 Lu 已把宏觀觀測、網絡傳播模型和 EnKF 型推斷結合,故本文不是第一次以聚合資料研究網絡疫情。他們的參數推斷問題,與這裏固定參數後比較分量包含率及寬度並不相同。Evensen 提供 ensemble 協方差方法的基礎,Niazi 與 Johansson 的疫情觀測器研究則補充模型相依狀態估計的脈絡。宏觀網絡同化EnKF 理論表述疫情狀態觀測器網絡 SIR 觀測器

所以,本文的價值是受控合成比較,以及解釋輸出究竟代表甚麼,不是提出新估計器、完整證明文獻優先權,或驗證真實疾病預測。下一步較有用的研究,應單獨隔離一個機制:在集合表示中保留更多跨區關係、比較獨立調校的機率濾波器,或改變資料何時及在哪裏可用。每個方向都需要自己的固定比較,不能只是為了改善本次成績而臨時調整。

結論

聚合報告可以限制隱藏疫情,但限制的意思取決於模型、量測日期及不確定性的表示方式。分布導向算法可以借助假設的協方差,抽取有用的地區結構;經驗證的集合算法可以保留所有允許軌跡,卻同時保留太多替代情況,以致不能支持精確預測。

閱讀這類圖像時,可以依次問三個問題:資料量度了哪個數量、哪一天?哪些假設讓圖上的區間有根據?把假設說清楚後,區間是否窄得足以回答科學問題?跳過最後一問,容易把安全的包絡誤當成有用答案;跳過前兩問,則容易把視覺上的信心誤當成知識。

本文不是公共衞生建議。模型假設固定速率、已知流動、已知延遲及有界合成誤差,沒有加入聲明以外的人口或通報機制。其教訓是方法論上的:不確定性分析既要保留方法漏掉真值的情況,也要保留它幾乎透過「甚麼都不排除」來包含真值的情況。

參考文獻

  1. Aronna, M. S. and Bliman, P.-A. (2017 preprint, revised 2018). Interval observer for uncertain time-varying SIR-SI epidemiological model of vector-borne disease.
  2. Meslem, N., Hably, A., Wang, Z. and Raïssi, T. (2024). Set-Valued State Estimator with Sparse and Delayed Measurements for Uncertain Discrete-Time Linear Systems.
  3. Efimov, D., Fridman, E., Polyakov, A., Perruquetti, W. and Richard, J.-P. (2016). Linear interval observers under delayed measurements and delay-dependent positivity.
  4. Rego, B. S., Scott, J. K., Raimondo, D. M. and Raffo, G. V. (2021). Set-valued state estimation of nonlinear discrete-time systems with nonlinear invariants based on constrained zonotopes.
  5. Gouzé, J.-L., Rapaport, A. and Hadj-Sadok, M. Z. (2000). Interval observers for uncertain biological systems.
  6. Niazi, M. U. B. and Johansson, K. H. (2022). Observer Design for the State Estimation of Epidemic Processes.
  7. Niazi, M. U. B. and Johansson, K. H. (2023). Parameterization-Free Observer Design for Nonlinear Systems: Application to the State Estimation of Networked SIR Epidemics.
  8. Wang, Y. and Lu, W. (2025). Estimating contagion dynamics models on networks via data assimilation.
  9. Evensen, G. (2003). The ensemble Kalman filter: Theoretical formulation and practical implementation.
  10. Nedialkov, N. S., Jackson, K. R. and Corliss, G. F. (1999). Validated solutions of initial value problems for ordinary differential equations.