機(jī)場(chǎng)建設(shè)總進(jìn)度計(jì)劃延期分布及貝葉斯估計(jì)

賈廣社，宋明禮，吳陸鋒，張普偉

（同濟(jì)大學(xué)經(jīng)濟(jì)與管理學(xué)院，上海200092）

建設(shè)工程總進(jìn)度計(jì)劃編制和跟蹤管理是指，在大型建設(shè)項(xiàng)目前期，編制一個(gè)指導(dǎo)整個(gè)工程設(shè)計(jì)準(zhǔn)備、設(shè)計(jì)、采購(gòu)、施工、動(dòng)用前準(zhǔn)備等方面進(jìn)度的綱要型進(jìn)度計(jì)劃，并定期跟蹤計(jì)劃的完成情況，提出管理建議［1］。施工進(jìn)度計(jì)劃往往由施工單位編制；而總進(jìn)度規(guī)劃不僅需要考慮施工進(jìn)度，還需要從項(xiàng)目整體角度考慮其他工作進(jìn)程及進(jìn)度風(fēng)險(xiǎn)因素，是由建設(shè)單位或建設(shè)單位委托專門的總進(jìn)度計(jì)劃管理單位進(jìn)行編制和跟蹤管理，總進(jìn)度計(jì)劃跟蹤管理屬于項(xiàng)目總控的工作范疇［1-2］。在我國(guó)近20年的機(jī)場(chǎng)建設(shè)工程實(shí)踐中，總進(jìn)度計(jì)劃的編制與跟蹤管理顯示出重要作用，對(duì)國(guó)內(nèi)14個(gè)大型機(jī)場(chǎng)新建、改擴(kuò)建工程的進(jìn)度總控咨詢成果（2000—2018年）顯示，總進(jìn)度計(jì)劃管理能夠幫助建設(shè)方有效把控項(xiàng)目進(jìn)度，滿足項(xiàng)目整體交付的剛性進(jìn)度要求。因此，實(shí)踐中更多的機(jī)場(chǎng)建設(shè)項(xiàng)目正在引入總進(jìn)度計(jì)劃管理。

總進(jìn)度計(jì)劃需要符合“戰(zhàn)略性、總體性、宏觀性”的要求，其制訂有別于傳統(tǒng)的完全自下而上的進(jìn)度計(jì)劃制訂方式［1，3-4］。王廣斌［1］闡述了大型建設(shè)項(xiàng)目總進(jìn)度綱要的概念、內(nèi)容、編制方法和步驟；針對(duì)大型機(jī)場(chǎng)航站樓建設(shè)，賈廣社等［3］通過(guò)總進(jìn)度計(jì)劃實(shí)踐成果資料及訪談，確定航站樓總控工序，進(jìn)而利用模擬方法估算航站樓建設(shè)總工期；羅晟等［4］通過(guò)上海世博會(huì)臨時(shí)場(chǎng)館項(xiàng)目群案例，闡釋了臨時(shí)場(chǎng)館項(xiàng)目群總進(jìn)度計(jì)劃編制與控制工作的基本框架、實(shí)施要點(diǎn)和標(biāo)準(zhǔn)流程；祁神軍等［5］在項(xiàng)目總控理論基礎(chǔ)之上，介紹從總進(jìn)度綱要直至執(zhí)行性計(jì)劃的大型建設(shè)項(xiàng)目總進(jìn)度計(jì)劃體系，并嘗試通過(guò)引入學(xué)習(xí)曲線，實(shí)現(xiàn)計(jì)劃的動(dòng)態(tài)調(diào)整。

機(jī)場(chǎng)建設(shè)項(xiàng)目總進(jìn)度計(jì)劃管理實(shí)踐中，將復(fù)雜的機(jī)場(chǎng)建設(shè)工程劃分為多個(gè)子系統(tǒng)，梳理出每個(gè)子系統(tǒng)工程各月度需要完成的各項(xiàng)具體工作；在總進(jìn)度計(jì)劃跟蹤管理過(guò)程中，管理者發(fā)現(xiàn)當(dāng)月未能按計(jì)劃完成的具體工作及其原因，制訂應(yīng)對(duì)策略，預(yù)測(cè)下月進(jìn)度風(fēng)險(xiǎn)情況，調(diào)整資源配置?？梢?jiàn)，有效預(yù)測(cè)各子系統(tǒng)工程工作延期（即進(jìn)度違約）發(fā)生數(shù)量，有利于有針對(duì)性地提前進(jìn)行風(fēng)險(xiǎn)預(yù)警和資源配置，從而保障總進(jìn)度計(jì)劃按期完成；也能為后續(xù)項(xiàng)目總進(jìn)度計(jì)劃的制訂提供參考依據(jù)。隨著總進(jìn)度計(jì)劃管理在機(jī)場(chǎng)建設(shè)工程實(shí)踐當(dāng)中的應(yīng)用，進(jìn)度違約歷史數(shù)據(jù)得到積累，如何利用這些數(shù)據(jù)，嘗試建立統(tǒng)計(jì)模型反映進(jìn)度風(fēng)險(xiǎn)發(fā)生規(guī)律，正是本研究所關(guān)注的問(wèn)題，此類統(tǒng)計(jì)模型的建立能夠幫助管理者從整體角度、系統(tǒng)角度刻畫和理解機(jī)場(chǎng)建設(shè)總進(jìn)度計(jì)劃實(shí)施過(guò)程中的進(jìn)度違約風(fēng)險(xiǎn)。

離散分布是常見(jiàn)的用來(lái)刻畫風(fēng)險(xiǎn)或違約發(fā)生可能的概率模型，除了在工業(yè)制造、交通規(guī)劃、消費(fèi)市場(chǎng)、金融信貸、保險(xiǎn)等研究領(lǐng)域的廣泛應(yīng)用之外，以泊松分布模型、負(fù)二項(xiàng)分布模型為代表的離散選擇模型也已經(jīng)在工程管理領(lǐng)域得到了應(yīng)用。Anastasopoulos等［6］運(yùn)用泊松回歸和負(fù)二項(xiàng)回歸模型研究了高速公路建設(shè)項(xiàng)目索賠發(fā)生次數(shù)及其影響因素；Erskine等［7］利用負(fù)二項(xiàng)分布描述城市下水基礎(chǔ)設(shè)施的故障發(fā)生頻數(shù)，并結(jié)合貝葉斯理論建立可更新的預(yù)測(cè)模型；Love等［8］利用負(fù)二項(xiàng)分布擬合建設(shè)工程中工人傷亡頻數(shù)和質(zhì)量缺陷發(fā)生頻數(shù)，并將二者作為建設(shè)工程項(xiàng)目安全和質(zhì)量方面風(fēng)險(xiǎn)的被動(dòng)預(yù)測(cè)指標(biāo)（passive lead-indicator）。動(dòng)態(tài)持續(xù)監(jiān)管進(jìn)度風(fēng)險(xiǎn)對(duì)提升項(xiàng)目績(jī)效有重要作用［9］?，F(xiàn)階段尚缺乏對(duì)進(jìn)度違約（default）頻數(shù)的研究，而工程的進(jìn)度違約同安全、質(zhì)量、成本等各類違約相似，都受到復(fù)雜不確定因素的影響，因而也適合用離散分布刻畫其發(fā)生頻數(shù)。所以，本文嘗試結(jié)合總進(jìn)度計(jì)劃管理實(shí)踐，構(gòu)建離散分布，刻畫每一個(gè)子系統(tǒng)一定時(shí)期內(nèi)工作延期發(fā)生頻數(shù)，從項(xiàng)目總體和宏觀層面反映機(jī)場(chǎng)建設(shè)工程項(xiàng)目各子系統(tǒng)的延期風(fēng)險(xiǎn)，為總進(jìn)度計(jì)劃風(fēng)險(xiǎn)預(yù)測(cè)和管理后續(xù)研究奠定基礎(chǔ)。

1 模型構(gòu)造前提與假設(shè)

（1）重大建設(shè)工程項(xiàng)目有其唯一性和獨(dú)特性，但機(jī)場(chǎng)建設(shè)工程作為一類具有嚴(yán)格專業(yè)標(biāo)準(zhǔn)的建設(shè)工程，其進(jìn)度特征如關(guān)鍵線路、延期風(fēng)險(xiǎn)等方面有一定相似性［10］。

（2）目前國(guó)內(nèi)大型機(jī)場(chǎng)建設(shè)工程的組織管理模式，均是以指揮部為最高管理機(jī)構(gòu)，下設(shè)飛行部、航站樓及衛(wèi)星廳部、設(shè)備部、計(jì)財(cái)部等專項(xiàng)部門和職能部門來(lái)進(jìn)行具體施工管理；且其工作分解結(jié)構(gòu)（WBS）較為相似。相應(yīng)地，建設(shè)工程子系統(tǒng)的劃分方法均類似于本案例，往往按照航站區(qū)、飛行區(qū)、綜合配套區(qū)等大的功能版塊，繼續(xù)細(xì)分成相應(yīng)的子系統(tǒng)，各子系統(tǒng)又包括若干能夠獨(dú)立執(zhí)行、監(jiān)測(cè)的工作活動(dòng)。以簡(jiǎn)化的某機(jī)場(chǎng)建設(shè)工程為例，其子系統(tǒng)劃分示例如圖1所示

（3）本文假設(shè)在該組織管理模式下，當(dāng)給予一定的總工期計(jì)劃目標(biāo)之后，在現(xiàn)有管理和技術(shù)水平之下，由于受到復(fù)雜不確定因素的影響，每個(gè)跟蹤期內(nèi)（本文為每月）每個(gè)子系統(tǒng)內(nèi)發(fā)生延期工作活動(dòng)數(shù)量服從一定形式的概率分布。

2 貝葉斯統(tǒng)計(jì)推斷與馬爾科夫蒙特卡洛（MCMC）方法

在對(duì)模型進(jìn)行分布參數(shù)估計(jì)時(shí)，由于大型工程進(jìn)度違約數(shù)據(jù)存在案例數(shù)量少，違約數(shù)量相對(duì)工作數(shù)量較小的特點(diǎn)，采用傳統(tǒng)參數(shù)估計(jì)方法效果并不理想，因此允許結(jié)合先驗(yàn)知識(shí)并能夠?qū)崟r(shí)更新的貝葉斯方法是一種較理想的估計(jì)方法；而借助MCMC方法，可以解決貝葉斯方法中的計(jì)算問(wèn)題，實(shí)現(xiàn)分布參數(shù)的貝葉斯統(tǒng)計(jì)推斷。

2.1 貝葉斯統(tǒng)計(jì)推斷與蒙特卡洛模擬

貝葉斯分析中，一個(gè)統(tǒng)計(jì)模型包含：先驗(yàn)分布π(θ)和似然函數(shù)L(x|θ)，其中x=(x1，x2，… ，xn)為容量為n的樣本，θ=(θ1，θ2，…，θn)為參數(shù)，樣本和參數(shù)都是隨機(jī)的。貝葉斯統(tǒng)計(jì)分析即是綜合先驗(yàn)信息、似然（由總體和樣本信息給出），得出后驗(yàn)信息。以連續(xù)分布情況為例，根據(jù)貝葉斯定理，θ的后驗(yàn)分布的密度函數(shù)為

式中：π(θ)為θ的先驗(yàn)密度函數(shù)；L(x|θ)為似然函數(shù)，θ∈Θ。

在平方損失函數(shù)下，參數(shù)θ的貝葉斯估計(jì)為

貝葉斯統(tǒng)計(jì)推斷就是要計(jì)算參數(shù)θ的后驗(yàn)分布的各階矩g(θ)（如后驗(yàn)均值、后驗(yàn)方差等）的貝葉斯估計(jì)量，即

可見(jiàn)，要得到參數(shù)θ及其各階矩g(θ)的后驗(yàn)分布，需要求以上公式中的積分，但在復(fù)雜情形或參數(shù)維數(shù)較大時(shí)，直接積分或數(shù)值積分難以實(shí)現(xiàn)?；诮y(tǒng)計(jì)計(jì)算理論的抽樣方法為貝葉斯統(tǒng)計(jì)分析中這類復(fù)雜積分的計(jì)算提供了解決方法。可以證明，如果抽樣自后驗(yàn)分布π(θ|x)的后驗(yàn)樣本為獨(dú)立的，根據(jù)大數(shù)定律，樣本均值-g依概率收斂于后驗(yàn)均值E[g(θ)|x]（即g(θ^)），只要樣本容量足夠大，估計(jì)精度可以達(dá)到任意所需要的精度，這就是蒙特卡洛（Monte Carlo，MC）估計(jì)。

2.2 MCMC方法簡(jiǎn)介

蒙特卡洛隨機(jī)模擬方法能夠?qū)崿F(xiàn)貝葉斯分析中復(fù)雜積分的運(yùn)算。但當(dāng)后驗(yàn)分布為高維、復(fù)雜且不常見(jiàn)的分布時(shí)，后驗(yàn)獨(dú)立樣本的抽取難度較大。馬爾科夫蒙特卡洛（Markov chain Monte Carlo，MCMC）方法便是利用馬爾科夫鏈從目標(biāo)分布π(θ|x)中抽取隨機(jī)樣本，進(jìn)而利用蒙特卡羅方法估算積分的方法。作為一種抽樣算法，相對(duì)于MC方法，MCMC提高了抽樣能力和效率。利用馬爾科夫鏈可將無(wú)法直接抽樣的復(fù)雜抽樣問(wèn)題轉(zhuǎn)化為簡(jiǎn)單抽樣問(wèn)題。具體而言，構(gòu)造一個(gè)平穩(wěn)分布恰好為后驗(yàn)分布π(θ|x)的馬爾科夫鏈，當(dāng)馬爾科夫鏈?zhǔn)諗亢?，即該鏈上的樣本點(diǎn)具有了一些良好性質(zhì)，如馬氏性、不可約性、遍歷性、非周期性之后，便可取該鏈上的樣本點(diǎn)序列作為后驗(yàn)“樣本”?？梢宰C明，這些“樣本”雖然并非相互獨(dú)立，但是卻與從π(θ|x)中直接抽取的獨(dú)立樣本具有同樣作用，可以用來(lái)以抽樣計(jì)算方法估算參數(shù)θ及其各階矩g(θ)。

3 數(shù)據(jù)來(lái)源及模型構(gòu)建

本研究數(shù)據(jù)來(lái)自浦東機(jī)場(chǎng)三期擴(kuò)建工程總進(jìn)度計(jì)劃管理實(shí)踐?？傔M(jìn)度計(jì)劃編制方法是：從子系統(tǒng)層面入手，結(jié)合專家經(jīng)驗(yàn)和宏觀規(guī)劃相關(guān)要求，通過(guò)政府相關(guān)單位訪談、建設(shè)單位及其職能部門訪談、各子系統(tǒng)主要承包商訪談，綜合制訂合理的、細(xì)分至子系統(tǒng)內(nèi)各項(xiàng)工作活動(dòng)的總進(jìn)度計(jì)劃體系；在跟蹤控制過(guò)程中，研究人員全過(guò)程參與業(yè)主主持召開的各類工作會(huì)議并形成會(huì)議紀(jì)要，實(shí)時(shí)進(jìn)行現(xiàn)場(chǎng)查勘獲取工程進(jìn)度資料和數(shù)據(jù)，每月根據(jù)會(huì)議內(nèi)容、現(xiàn)場(chǎng)資料制訂總進(jìn)度月報(bào)，并進(jìn)行風(fēng)險(xiǎn)預(yù)警。所截取的進(jìn)度計(jì)劃跟蹤數(shù)據(jù)來(lái)自2016年7月至2018年7月，總計(jì)25個(gè)月。所選取的13個(gè)子系統(tǒng)（如表1所示）在該時(shí)間段內(nèi)均已正式展開系統(tǒng)內(nèi)建設(shè)工作。

各個(gè)子系統(tǒng)工程的進(jìn)度計(jì)劃跟蹤管理表的示例如表2所示（省略部分未使用信息），該月份計(jì)劃工作活動(dòng)57項(xiàng)，實(shí)際完成53項(xiàng)，發(fā)生延期4項(xiàng)。本文使用的數(shù)據(jù)包括月度各系統(tǒng)內(nèi)開展工作活動(dòng)的總數(shù)量，以及其中發(fā)生延期的工作數(shù)量。

表1 案例機(jī)場(chǎng)項(xiàng)目所選取的子系統(tǒng)Tab.1 Subsystems of the case airport construction project

對(duì)于每個(gè)子系統(tǒng)，構(gòu)建一個(gè)離散概率分布，來(lái)刻畫一定時(shí)期內(nèi)可能發(fā)生的延期活動(dòng)數(shù)量（延期活動(dòng)頻數(shù)）。令延期活動(dòng)頻數(shù)分布表示為f(nt|θ)，其中nt表示一個(gè)子系統(tǒng)在第t個(gè)月發(fā)生的延期活動(dòng)數(shù)量，由于所跟蹤的時(shí)間共m個(gè)月，所以t=1，2，…，m。θ表示該分布的參數(shù)，則某一子系統(tǒng)在m個(gè)月內(nèi)發(fā)生延期活動(dòng)數(shù)量的似然函數(shù)為

對(duì)于一定時(shí)間內(nèi)超期活動(dòng)發(fā)生若干次數(shù)的概率，適宜采用描述離散隨機(jī)事件的概率分布來(lái)進(jìn)行刻畫，因此本文中選用泊松分布（Poisson distribution）和 負(fù) 二 項(xiàng) 分 布（negative binomial distribution）進(jìn)行研究。

4 參數(shù)估計(jì)過(guò)程

4.1 極大似然方法

對(duì)于泊松分布和負(fù)二項(xiàng)分布假設(shè)，根據(jù)案例數(shù)據(jù)，采用極大似然估計(jì)來(lái)計(jì)算分布參數(shù)，具體采用牛頓法迭代求解似然函數(shù)。表3展示了延期發(fā)生數(shù)量概率分布參數(shù)的極大似然估計(jì)值。表中λ為泊松分布參數(shù)，θ、r為負(fù)二項(xiàng)分布參數(shù)。

4.2 貝葉斯方法

4.2.1 泊松分布假設(shè)

對(duì)于每一個(gè)子系統(tǒng)，假設(shè)其在25個(gè)月中，每個(gè)月發(fā)生的延期工作數(shù)量nt服從泊松分布，即

式中：t=1，2，…，25。

則似然函數(shù)形式為

在貝葉斯后驗(yàn)分析當(dāng)中，先驗(yàn)分布通常選取共軛先驗(yàn)分布形式，這是因?yàn)樨惾~斯后驗(yàn)分布既反映過(guò)往的經(jīng)驗(yàn)（參數(shù)的先驗(yàn)分布），又反映樣本的信息（觀測(cè)數(shù)據(jù)），而共軛型分布要求先驗(yàn)分布與后驗(yàn)分布為同一種分布類型，所以選取共軛先驗(yàn)分布就意味著使經(jīng)驗(yàn)知識(shí)與現(xiàn)有樣本信息具有某種同一性，它們能夠轉(zhuǎn)化為同一類的經(jīng)驗(yàn)知識(shí)；并且，所得到后驗(yàn)分布作為新的經(jīng)驗(yàn)知識(shí)進(jìn)行下一步貝葉斯信念更新后，所得到的新的后驗(yàn)分布仍然是同一類型的。對(duì)于泊松分布的參數(shù)λ，取其共軛先驗(yàn)分布Gamma分布作為先驗(yàn)分布，即

式中：a，b為超參數(shù)。

超參數(shù)可根據(jù)先驗(yàn)信息確定，本研究假設(shè)總進(jìn)度計(jì)劃實(shí)施受到難以預(yù)知的復(fù)雜隨機(jī)因素影響，采用提供十分有限信息的非信息（uninformative）先驗(yàn)。其一，根據(jù)課題組以往參加的機(jī)場(chǎng)建設(shè)工程實(shí)踐記錄，在類似項(xiàng)目分類以及工作實(shí)施過(guò)程中，每月發(fā)生的延期工作數(shù)量不會(huì)超過(guò)10個(gè)，除此之外，缺乏其他先驗(yàn)信息，由此可采用先驗(yàn)分布λ～Γ(5，1)；其二，采用模糊先驗(yàn)（后文將進(jìn)行論述）。需要指出的是，本研究假設(shè)所選取的先驗(yàn)分布，是提供了有限的、較為“模糊”的客觀信息。這兩類先驗(yàn)都可稱為非信息先驗(yàn)（uninformative prior），這是相對(duì)于信息先驗(yàn)（informative prior）而言，非信息先驗(yàn)表示對(duì)于研究對(duì)象的先驗(yàn)知識(shí)十分有限。而信息先驗(yàn)則是指通過(guò)實(shí)踐經(jīng)驗(yàn)、專家知識(shí)等形成了對(duì)于研究對(duì)象較為豐富的經(jīng)驗(yàn)知識(shí)。本文采用非信息先驗(yàn)，旨在介紹貝葉斯統(tǒng)計(jì)推斷方法在延期頻數(shù)分布研究中的應(yīng)用，并給出具體實(shí)現(xiàn)過(guò)程示例；而至于先驗(yàn)信息的進(jìn)一步收集、豐富和引入，則可按相應(yīng)的理論和方法進(jìn)行進(jìn)一步研究。

表2 某月“航站樓/衛(wèi)星廳工程”子系統(tǒng)內(nèi)計(jì)劃活動(dòng)及完成情況Tab.2 Planned activities and completion of subsystem“terminal/satellite hall project”in a certain month

表3 頻數(shù)分布參數(shù)的極大似然估計(jì)值Tab.3 Maximum likelihood estimates of frequency distributions

先驗(yàn)分布λ～Γ(5,1)的圖像如圖2所示。

圖2 參數(shù)λ的共軛先驗(yàn)分布Γ（5，1）Fig.2 Prior Γ（5，1）for parameter λ

由于所選取的是共軛先驗(yàn)分布，參數(shù)λ的后驗(yàn)分布也將服從Gamma分布形式，因此可以直接采用MCMC算法中的Gibbs抽樣方法進(jìn)行抽樣計(jì)算，得到后驗(yàn)樣本，進(jìn)行后驗(yàn)分析。以3號(hào)子系統(tǒng)“捷運(yùn)系統(tǒng)工程”為例圖3為3號(hào)子系統(tǒng)泊松分布參數(shù)λ的后驗(yàn)抽樣樣本路徑。可以看出，圖3中無(wú)持續(xù)扁平區(qū)域，且波動(dòng)并不過(guò)分集中。

圖3 泊松分布參數(shù)的后驗(yàn)抽樣樣本路徑（3號(hào)子系統(tǒng)）Fig.3 Parameter sample paths for Poisson distribution(Subsystem 3)

類似地，從圖4所示的后驗(yàn)抽樣跡圖可見(jiàn)，MCMC算法模擬產(chǎn)生的參數(shù)值無(wú)明顯趨勢(shì)特征（tendency）或周期特征（periodicity），所產(chǎn)生樣本參數(shù)值絕大部分都在同一值域之內(nèi)。說(shuō)明抽樣確實(shí)符合馬爾科夫鏈特征，即每一狀態(tài)都與其歷史狀態(tài)相關(guān)；同時(shí)每個(gè)值隨機(jī)持續(xù)波動(dòng)，符合蒙特卡洛方法獨(dú)立抽樣的需求。

圖4 泊松分布參數(shù)的后驗(yàn)抽樣跡圖（3號(hào)子系統(tǒng)）Fig.4 Parameter sample trace plot for Poisson distribution(Subsystem 3)

以上兩步為從不同角度進(jìn)行的算法收斂性檢驗(yàn)，除此之外，還可以用遍歷均值圖（ergodic mean）方法對(duì)收斂性進(jìn)行檢驗(yàn)。從非周期性、馬氏性、遍歷性等不同角度來(lái)判斷，此處所使用的MCMC算法皆收斂良好。針對(duì)本文選取的參數(shù)先驗(yàn)分布及其超參數(shù)，該抽樣算法能夠成功地在參數(shù)取值空間隨機(jī)抽取樣本［11］，且并沒(méi)有長(zhǎng)時(shí)間停滯在參數(shù)空間的某一個(gè)小區(qū)域進(jìn)行抽樣），因此所產(chǎn)生的馬爾科夫鏈上的抽樣點(diǎn)雖然并非相互獨(dú)立，但是卻與從后驗(yàn)分布中直接抽取的獨(dú)立樣本具有同樣作用，可以用來(lái)以抽樣計(jì)算方法估算參數(shù)θ。注意，這里舍棄了前20 000個(gè)樣本點(diǎn)，作為抽樣預(yù)熱階段（burn-in period），保留下來(lái)的樣本點(diǎn)是更加接近馬爾科夫鏈平穩(wěn)分布的樣本點(diǎn)，本文據(jù)此進(jìn)行后驗(yàn)分布統(tǒng)計(jì)量的計(jì)算。泊松分布參數(shù)λ的后驗(yàn)密度估計(jì)如圖5所示。

對(duì)于3號(hào)子系統(tǒng)“捷運(yùn)系統(tǒng)工程”，泊松分布假設(shè)下MCMC過(guò)程產(chǎn)生參數(shù)樣本的自相關(guān)圖如圖6所示，其中Lag值為滯后階數(shù)，如Lag值為3時(shí)，表示該抽樣樣本與其后第3個(gè)抽樣樣本之間的自相關(guān)程度。表4則展示了在Lag值取1至49時(shí)，參數(shù)λ抽樣結(jié)果具體的自相關(guān)系數(shù)值。從圖6與表4可以看出，后驗(yàn)抽樣樣本的自相關(guān)程度非常低，從判斷抽樣算法收斂的角度來(lái)看，這是一個(gè)理想結(jié)果。

圖5 泊松分布參數(shù)的后驗(yàn)密度估計(jì)圖（3號(hào)子系統(tǒng)，后驗(yàn)抽樣樣本數(shù)為80 000）Fig.5 Parameter density plot for Poisson distribution(Subsystem 3,sample:80 000)

圖6 泊松分布MCMC過(guò)程的自相關(guān)圖（3號(hào)子系統(tǒng)）Fig.6 Autocorrelation plot for Poisson distribution(Subsystem 3)

蒙特卡洛誤差（MC error）Me是基于MCMC方法輸出結(jié)果所計(jì)算的重要指標(biāo)，用于衡量模擬計(jì)算而得的參數(shù)值的變動(dòng)程度。本文采用分塊均值（batch mean）計(jì)算Me值，簡(jiǎn)而言之，就是將MCMC方法產(chǎn)生的所有樣本T均分為K個(gè)分塊，每個(gè)分塊含有v=T/K個(gè)樣本，先計(jì)算每個(gè)分塊的均值，對(duì)所得的所有分塊均值求標(biāo)準(zhǔn)差，即為分塊均值法計(jì)算的Me值?？梢?jiàn)，較小的Me值意味著根據(jù)模擬樣本所計(jì)算的參數(shù)估計(jì)值變動(dòng)較小，即對(duì)參數(shù)的估計(jì)具有較高的精確度［11］。并且，Me值也是監(jiān)測(cè)MCMC算法收斂性地重要指標(biāo)，較小的Me值是算法收斂的必要條件。

各子系統(tǒng)后驗(yàn)抽樣過(guò)程所計(jì)算的Me值如表5所示。從表5可見(jiàn)，相對(duì)于后驗(yàn)樣本標(biāo)準(zhǔn)差，Me值均很小，可以說(shuō)明對(duì)參數(shù)的估計(jì)比較準(zhǔn)確。

本案例中，對(duì)其他各個(gè)子系統(tǒng)參數(shù)進(jìn)行抽樣后，以上各后驗(yàn)分析過(guò)程也都能得出類似結(jié)果，證明MCMC算法收斂，且估計(jì)結(jié)果符合精度要求。

另外，如前所述，還可采用“模糊先驗(yàn)”［12］（vague prior）作為先驗(yàn)分布，具體來(lái)說(shuō)，令Gamma分布Γ(a，b)的均值和方差分別等于先驗(yàn)均值-λ和先驗(yàn)方差S2λ，由于Gamma分布Γ(a，b)的均值和方差為

表4 泊松分布MCMC過(guò)程的自相關(guān)值（3號(hào)子系統(tǒng)）Tab.4 Autocorrelations for Poisson distribution(Subsystem 3)

表5 各子系統(tǒng)泊松分布MCMC過(guò)程的MC error值Tab.5 MC errors for Poisson distribution for subsystems

則其參數(shù)a，b可由下式給出：

4.2.2 負(fù)二項(xiàng)分布假設(shè)

對(duì)于每一個(gè)子系統(tǒng)，假設(shè)其在25個(gè)月中，每個(gè)月發(fā)生的延期工作數(shù)量nt服從負(fù)二項(xiàng)分布，即

式中：t=1，2，…，25。

則似然函數(shù)形式為

根據(jù)負(fù)二項(xiàng)分布的定義，參數(shù)θ為伯努利試驗(yàn)“成功”的概率，而r則是表示在該負(fù)二項(xiàng)分布中“成功”的試驗(yàn)數(shù)量；對(duì)于負(fù)二項(xiàng)分布參數(shù)θ，選取的共軛先驗(yàn)分布為Beta分布β(a，b)。關(guān)于其超參數(shù)a，b的確定，亦采用以下兩種方式。

首先，使用類似于前述模糊先驗(yàn)，即以往研究中多采用的先驗(yàn)矩計(jì)算超參數(shù)a，b。具體來(lái)說(shuō)，令Beta分布β(a，b)的均值和方差分別等于先驗(yàn)均值和先驗(yàn)方差，即

則可解得其超參數(shù)a，b的估計(jì)分別為

另外，還可采用無(wú)信息先驗(yàn)分布確定超參數(shù)a，b。無(wú)信息先驗(yàn)分布是指除了參數(shù)的取值范圍和參數(shù)在總體分布中的地位之外，再也不包含參數(shù)的任何信息的先驗(yàn)分布?？梢宰C明，Jefferys準(zhǔn)則下，分布參數(shù)θ的無(wú)信息先驗(yàn)分布為Beta分布β(0.5，0.5)；而Lindley準(zhǔn)則下，參數(shù)θ的無(wú)信息先驗(yàn)分布為Beta分布β(0，0)。事實(shí)上，從貝葉斯估計(jì)風(fēng)險(xiǎn)角度看，β(1，1)，β(0.5，0.5)及β(0，0)三種無(wú)信息先驗(yàn)分布都具有較小估計(jì)風(fēng)險(xiǎn)［13］，都是合理可用的。目前無(wú)論是統(tǒng)計(jì)理論研究還是應(yīng)用研究，無(wú)信息先驗(yàn)分布的采用越來(lái)越多，被認(rèn)為是“客觀”和可接受的，這也是貝葉斯統(tǒng)計(jì)學(xué)派近幾十年研究中最成功的部分之一［14］。在這里，由于Beta分布刻畫的是一個(gè)取值在（0，1）的參數(shù)，則其方差不宜過(guò)大，同時(shí)，固定超參數(shù)a，b的取值，構(gòu)成U型Beta先驗(yàn)分布，以此表示先驗(yàn)信息的有限［12］，因此本文選用Beta(0.5，0.5)作為負(fù)二項(xiàng)分布參數(shù)θ的無(wú)信息先驗(yàn)，即

對(duì)于負(fù)二項(xiàng)分布的參數(shù)r，所選取的先驗(yàn)分布為泊松分布，并且根據(jù)負(fù)二項(xiàng)分布定義，參數(shù)r應(yīng)當(dāng)為正整數(shù)，因此，可令

其中t(v)表示對(duì)于v取整。采用前述模糊先驗(yàn)，本研究中利用案例數(shù)據(jù)的先驗(yàn)矩計(jì)算超參數(shù)h。具體來(lái)說(shuō)，令泊松分布P(h)的均值等于先驗(yàn)均值-r，即可求得超參數(shù)h。最終確定的超參數(shù)h取值如表6所示。

表6 各子系統(tǒng)負(fù)二項(xiàng)分布先驗(yàn)下的超參數(shù)h取值Tab.6 Hyperparameterhfor negative binomial prior

以4號(hào)子系統(tǒng)“衛(wèi)星廳港灣機(jī)坪飛行區(qū)工程”為例，經(jīng)MCMC抽樣算法（具體為Metropolis Hastings算法）抽樣之后，參數(shù)θ與r樣本路徑如圖7所示。

參數(shù)θ與r的后驗(yàn)抽樣樣本路徑較為特殊的原因在于，參數(shù)r必須取整數(shù)值，而在抽樣算法中參數(shù)θ與r是聯(lián)動(dòng)變化，同時(shí)進(jìn)行更新迭代的，因此其路徑形成如圖7所示形態(tài)。

圖7 負(fù)二項(xiàng)分布參數(shù)的后驗(yàn)抽樣樣本路徑（4號(hào)子系統(tǒng)）Fig.7 Parameter sample paths for negative binomial distribution(Subsystem 4)

圖8 和圖9分別為去掉前20 000個(gè)burn-in過(guò)程樣本后，參數(shù)θ的邊際后驗(yàn)密度圖與參數(shù)r的邊際后驗(yàn)直方圖。一般情況下，多峰的密度圖可能是由于馬爾科夫鏈未收斂，或者與數(shù)據(jù)相沖突的過(guò)強(qiáng)先驗(yàn)假設(shè)所導(dǎo)致［11］，然而在本例中，多峰現(xiàn)象則是由于負(fù)二項(xiàng)分布參數(shù)θ與r的內(nèi)在聯(lián)系，因此算法當(dāng)中θ與r聯(lián)動(dòng)構(gòu)造，同時(shí)進(jìn)行更新迭代?？梢钥吹?，參數(shù)θ的后驗(yàn)密度圖的兩個(gè)波峰（modal value），分別與參數(shù)r的后驗(yàn)頻率分布直方圖中頻率最高的兩個(gè)r值相對(duì)應(yīng)。圖10為參數(shù)θ與r的自相關(guān)圖（ACF圖），兩者皆呈現(xiàn)出一個(gè)隨著Lag值增大，相對(duì)緩慢下降，最終趨近于0的過(guò)程。相似的下降過(guò)程也反映出Metropolis Hastings抽樣算法構(gòu)造當(dāng)中參數(shù)θ與r的關(guān)聯(lián)，趨近于0的自相關(guān)函數(shù)值是算法收斂的必要條件。

圖8 負(fù)二項(xiàng)分布參數(shù)θ的后驗(yàn)密度估計(jì)圖（4號(hào)子系統(tǒng)）Fig.8Density plot for parameter θ (Subsystem 4)

圖9 負(fù)二項(xiàng)分布參數(shù)r的后驗(yàn)頻率分布直方圖（4號(hào)子系統(tǒng)）Fig.9 Frequency distribution histogram for parameter r(Subsystem 4)

基于本案例數(shù)據(jù)，對(duì)于各個(gè)子系統(tǒng)，參數(shù)θ與r后驗(yàn)抽樣過(guò)程中計(jì)算得出的MC errors以及遍歷均值（ergodic mean）結(jié)果，能夠共同證明算法的收斂性，分析過(guò)程類似上文所述，這里不再展示。最終，去掉前20 000個(gè)burn-in過(guò)程樣本后，以80 000個(gè)模擬抽樣樣本計(jì)算，分別得到了在泊松分布假設(shè)與負(fù)二項(xiàng)分布假設(shè)之下，參數(shù)的貝葉斯后驗(yàn)估計(jì)值（后驗(yàn)均值）如表7所示。

圖10 負(fù)二項(xiàng)分布MCMC過(guò)程自相關(guān)圖（4號(hào)子系統(tǒng)）Fig.10 Autocorrelation plot for negative binomial distribution(Subsystem 4)

表7 各子系統(tǒng)延期頻數(shù)分布參數(shù)的貝葉斯估計(jì)值Tab.7 Parameters of frequency distributions estimated with Bayesian approach

5 結(jié)果分析

以上分布參數(shù)貝葉斯估計(jì)詳細(xì)過(guò)程，給出了先驗(yàn)類型選擇、超參數(shù)設(shè)定、算法收斂判斷等貝葉斯估計(jì)步驟及結(jié)果，證明了其可行性。由圖11可知，從貝葉斯參數(shù)估計(jì)結(jié)果和極大似然（MLE）參數(shù)估計(jì)結(jié)果對(duì)比來(lái)看，在當(dāng)前數(shù)據(jù)下，各子系統(tǒng)間的MLE估計(jì)結(jié)果差別較大，而貝葉斯方法減少了估計(jì)結(jié)果的波動(dòng)性，這與金融領(lǐng)域類似應(yīng)用研究［15］的結(jié)論相一致，即相對(duì)于極大似然估計(jì)，結(jié)合歷史數(shù)據(jù)和先驗(yàn)信息的貝葉斯估計(jì)結(jié)果在低違約率事件違約頻數(shù)建模和預(yù)測(cè)方面具有穩(wěn)健性優(yōu)勢(shì)。機(jī)場(chǎng)建設(shè)工程總進(jìn)度數(shù)據(jù)也具有獲取難度大，案例數(shù)量有限的特征，因此在使用離散分布刻畫進(jìn)度違約可能時(shí)，分布參數(shù)的貝葉斯估計(jì)方法應(yīng)當(dāng)作為傳統(tǒng)估計(jì)方法的必要補(bǔ)充。這有助于管理者更好地理解機(jī)場(chǎng)建設(shè)總進(jìn)度計(jì)劃進(jìn)度違約發(fā)生規(guī)律。

圖11 延期發(fā)生頻數(shù)分布參數(shù)MLE估計(jì)值與貝葉斯估計(jì)值Fig.11 MLE and MCMC estimates for frequency distributions

采用泊松分布假設(shè)估計(jì)分布參數(shù)時(shí)，貝葉斯估計(jì)過(guò)程采用了兩種先驗(yàn)，得到的結(jié)果十分相似（圖11a），但是其后驗(yàn)估計(jì)結(jié)果都與MLE結(jié)果不相類似。而相比之下，負(fù)二項(xiàng)分布參數(shù)的貝葉斯估計(jì)結(jié)果與MLE結(jié)果比較接近，這表現(xiàn)在各個(gè)子系統(tǒng)參數(shù)的絕對(duì)值和相對(duì)值上（圖11b、11c）。這說(shuō)明在應(yīng)用泊松分布假設(shè)時(shí)還需要更多準(zhǔn)確先驗(yàn)知識(shí)的使用，完善泊松分布形式、參數(shù)及其先驗(yàn)的設(shè)定，例如采用其他形式的超先驗(yàn)分布，或新加入超參數(shù)去納入更多的先驗(yàn)信息［11］。相比較而言，利用負(fù)二項(xiàng)分布先驗(yàn)，貝葉斯估計(jì)結(jié)果與MLE結(jié)果較為接近。若假設(shè)該工程樣本數(shù)據(jù)具有代表性，則說(shuō)明負(fù)二項(xiàng)分布能夠較好地刻畫超期工作發(fā)生頻數(shù)，觀察進(jìn)度違約月度數(shù)據(jù)，其方差相對(duì)均值較大，存在過(guò)離散（overdispersion）趨勢(shì)，這也是負(fù)二項(xiàng)分布能夠較好適用的原因之一［7，15］。另外，本案例中觀察到的樣本數(shù)據(jù)可能只部分地反應(yīng)機(jī)場(chǎng)建設(shè)工程總進(jìn)度計(jì)劃延期工作發(fā)生的客觀規(guī)律，那么需要更多后續(xù)案例的數(shù)據(jù)來(lái)進(jìn)行建模和驗(yàn)證（validation）［16］，而可更新性也正是貝葉斯模型的優(yōu)勢(shì)之一。

6 結(jié)語(yǔ)

基于工程實(shí)踐提出一種刻畫機(jī)場(chǎng)建設(shè)工程總進(jìn)度計(jì)劃延期發(fā)生風(fēng)險(xiǎn)的模型；并利用實(shí)踐歷史數(shù)據(jù)，通過(guò)MCMC方法的實(shí)施說(shuō)明貝葉斯分布參數(shù)估計(jì)的可行性和操作過(guò)程，估計(jì)結(jié)果證明了貝葉斯方法的穩(wěn)健性優(yōu)勢(shì)和必要性。該分布模型可以幫助管理者理解在現(xiàn)有機(jī)場(chǎng)建設(shè)項(xiàng)目總進(jìn)度計(jì)劃工作分解框架及工期目標(biāo)之下，不同子系統(tǒng)中的工作活動(dòng)發(fā)生延期的可能性，并可根據(jù)實(shí)際發(fā)生數(shù)據(jù)不斷更新模型。模型有助于管理者對(duì)按期完成難度較大的工作適度增加關(guān)注和投入，保證總進(jìn)度計(jì)劃順利實(shí)施；在新制訂總進(jìn)度計(jì)劃時(shí)，關(guān)于延期發(fā)生可能的信息也有助于給各子系統(tǒng)分配合理的工期。例如，在目標(biāo)允許的情況下，給延期發(fā)生可能較高的子系統(tǒng)適當(dāng)增加工期，提高計(jì)劃按時(shí)完成的概率，避免延期引起的超額投入和其他問(wèn)題。需要指出的是，本研究所構(gòu)建的基于風(fēng)險(xiǎn)隨機(jī)假設(shè)、根據(jù)工程實(shí)踐定量數(shù)據(jù)、宏觀刻畫進(jìn)度風(fēng)險(xiǎn)可能性的概率模型，應(yīng)當(dāng)作為機(jī)場(chǎng)建設(shè)項(xiàng)目總進(jìn)度計(jì)劃制訂、管理和風(fēng)險(xiǎn)預(yù)測(cè)的輔助手段。但與概率風(fēng)險(xiǎn)模型在其他任何行業(yè)的應(yīng)用相同，其在機(jī)場(chǎng)建設(shè)項(xiàng)目中的應(yīng)用也必須結(jié)合定性研究結(jié)論和其他手段，以保證有效性。例如進(jìn)度違約的具體原因、子系統(tǒng)間的相互制約、延期工作的改進(jìn)方向等微觀具體問(wèn)題需要進(jìn)一步案例研究和實(shí)踐總結(jié)來(lái)揭示。后續(xù)研究中，通過(guò)關(guān)于先驗(yàn)知識(shí)獲取及先驗(yàn)分布確定的進(jìn)一步研究，可以將更多定性數(shù)據(jù)、專家知識(shí)引入模型［12］；通過(guò)建立泊松回歸、負(fù)二項(xiàng)回歸等回歸模型，可實(shí)現(xiàn)大型建設(shè)工程各總進(jìn)度計(jì)劃超期風(fēng)險(xiǎn)影響因素的判別；另外，在掌握超期頻數(shù)的基礎(chǔ)上，結(jié)合超期所造成損失的嚴(yán)重程度，可以實(shí)現(xiàn)對(duì)進(jìn)度風(fēng)險(xiǎn)更進(jìn)一步的刻畫和分析［17-19］。