基于灰色權(quán)重馬爾可夫鏈工傷事故預(yù)測(cè)

張昌爽

（重慶大學(xué)資源與安全學(xué)院，重慶 400044）

工傷事故灰色預(yù)測(cè)模型展現(xiàn)工傷事故時(shí)間狀態(tài)序列數(shù)據(jù)及其相應(yīng)變化總體波動(dòng)趨勢(shì)。建立灰色馬爾可夫預(yù)測(cè)模型，用馬爾可夫預(yù)測(cè)法確定未來(lái)狀態(tài)轉(zhuǎn)移規(guī)律，馬爾可夫灰色預(yù)測(cè)具備顯著優(yōu)勢(shì)。馬爾可夫鏈為狀態(tài)可列、時(shí)間齊次馬爾可夫過(guò)程[1-3]。

1 灰色權(quán)重馬爾可夫鏈預(yù)測(cè)步驟

1.1 狀態(tài)劃分

由有序聚類法構(gòu)造工傷人數(shù)分級(jí)標(biāo)準(zhǔn)，構(gòu)建馬爾可夫鏈狀態(tài)指標(biāo)值分級(jí)空間，確定各時(shí)段指標(biāo)所對(duì)應(yīng)狀態(tài)[4-6]。變量x1，x2，x3歸類為{xi,…,xj},j≥i，定義均值向量：

式中，D(i,j)為量化變量段內(nèi)部各變量間的差異值，以徑向值{xi，…，xj}(j≥i)進(jìn)行確定；徑向量數(shù)值越小，表示區(qū)間向量差異越小，彼此越接近。

將n個(gè)變量序列分為K類，分類誤差函數(shù)如下：

式中，k未知分類數(shù)，計(jì)算k=1，2，…，n所有最佳分割。

誤差函數(shù)e[P(n,k)]通過(guò)線性規(guī)劃，對(duì)應(yīng)最優(yōu)分類標(biāo)準(zhǔn)通過(guò)在min{e[P(n,k)]}中得以確定，確定分類數(shù)K變化趨勢(shì)，作出誤差函數(shù)e[P(n,k)]與K的關(guān)系圖，函數(shù)曲線拐點(diǎn)對(duì)應(yīng)的K值為最優(yōu)分類數(shù)。

1.2 確定狀態(tài)轉(zhuǎn)移概率矩陣

轉(zhuǎn)移概率以{x1，x2，…，xn}定義馬爾可夫鏈指標(biāo)序列，包含m個(gè)狀態(tài)，指標(biāo)值序列從狀態(tài)i經(jīng)過(guò)1步轉(zhuǎn)移到達(dá)狀態(tài)j的頻數(shù)以fij表示，i，j∈E，轉(zhuǎn)移步長(zhǎng)是1個(gè)時(shí)間單位，2個(gè)甚至m個(gè)時(shí)間單位。

式中，fij為轉(zhuǎn)移頻數(shù)矩陣(fij)組成單元，轉(zhuǎn)移頻數(shù)矩陣的第i行第j列元素fij除以對(duì)應(yīng)各行的單元值總和所得，表示為Pij，i，j∈E。

1.3 馬氏性檢驗(yàn)

用χ2統(tǒng)計(jì)量檢驗(yàn)離散型序列隨機(jī)變量。M表示隨機(jī)變量指標(biāo)值序列存在的演變狀態(tài)，fij為指標(biāo)值序列{x1，x2，…，xn}從狀態(tài)i經(jīng)過(guò)1步轉(zhuǎn)移到達(dá)狀態(tài)j的頻數(shù)值，i，j∈E。

式中，p.j為邊際概率，轉(zhuǎn)移頻數(shù)矩陣第j列之和除以矩陣所有行和列總和。

當(dāng)n的基數(shù)足夠大，統(tǒng)計(jì)量驗(yàn)證指標(biāo)：

式中，pij為轉(zhuǎn)移概率。

在顯著性水平a，服從自由度(m-1)2的χ2分布，查表得分位點(diǎn)的客觀數(shù)值，計(jì)算后得研究指標(biāo)統(tǒng)計(jì)量χ2值。若，則xi符合馬氏性檢驗(yàn)標(biāo)準(zhǔn)，否則該數(shù)據(jù)序列不適用于馬爾可夫鏈來(lái)預(yù)測(cè)。

1.4 各階自相關(guān)系數(shù)計(jì)算

自相關(guān)系數(shù)rk，k∈e；rk為第k階（滯時(shí)為k季度）自相關(guān)系數(shù)；xl為第l 季度工傷人數(shù)；x季度平均受傷人數(shù)均值；n季度平均受傷人數(shù)序列長(zhǎng)度。

式中，rk為各階自相關(guān)系數(shù)。

規(guī)范各階自相關(guān)系數(shù)，即：

式中，wk步長(zhǎng)馬爾可夫鏈權(quán)重；m預(yù)測(cè)最大階數(shù)。

以狀態(tài)特征值法、線性插值法予以推斷，將預(yù)測(cè)數(shù)據(jù)加入原矩陣序列之中，重復(fù)上述步驟。

式中，μ為狀態(tài)特征值。

其中，i≤μ＜i+1，β＞0為狀態(tài)調(diào)整因子，狀態(tài)區(qū)間下限為，上限為，預(yù)測(cè)指標(biāo)值通過(guò)上下限兩類算法，結(jié)果一致。

式中，x為指標(biāo)預(yù)測(cè)值；k=0，狀態(tài)值增加與指標(biāo)值一致；否則，k=1。

非零向量πj={π1,π2,…,πE}，使得πjPij=πj，其中Pij為狀態(tài)轉(zhuǎn)移概率矩陣，有πj={π1,π2,…,πE}為馬爾可夫鏈平穩(wěn)分布，處于平穩(wěn)狀態(tài)。當(dāng)存在m＞0，pm中諸元素皆非負(fù)非零，則p為正規(guī)概率矩陣。

2 模型應(yīng)用

以2010年1季度至2019年第二季度共38季度工傷人數(shù)預(yù)測(cè)2019年3季度的受傷人數(shù)，將2019年3季度實(shí)際預(yù)測(cè)數(shù)據(jù)加入序列模型中，預(yù)測(cè)2019年4季度受傷人數(shù)，并與社保局實(shí)際備案的工傷人數(shù)對(duì)比，驗(yàn)證預(yù)測(cè)效果。

2.1 分級(jí)準(zhǔn)則確定狀態(tài)

用有序聚類法將序列分5個(gè)人數(shù)區(qū)間比較，工傷人數(shù)等級(jí)標(biāo)準(zhǔn)以1～5表示，狀態(tài)依次由低到高（見(jiàn)表1）。

表1 工傷人數(shù)狀態(tài)劃分表

2.2 馬氏性檢驗(yàn)

2010年1 季度至2019年第4季度共40個(gè)季度某企業(yè)工傷事故人數(shù)序列作馬氏性檢驗(yàn)。每季度工傷人數(shù)頻數(shù)轉(zhuǎn)移矩陣(fij)5×5、步長(zhǎng)為1的1步轉(zhuǎn)移概率矩陣(pij)5×5如下：

由(fij)5×5及(pij)5×5得邊際概率、統(tǒng)計(jì)量χ2值（見(jiàn)表2、表3）：

表2 邊際概率

表3 統(tǒng)計(jì)量χ2計(jì)算表

計(jì)算χ2為54.5556，在顯著性水平a=0.05下，查表得分位點(diǎn)由于因此工傷人數(shù)序列滿足馬氏性檢驗(yàn)。

2.3 馬爾可夫鏈權(quán)重確定

2010 年1 季度至2019 年第二季度共38 季度工傷數(shù)據(jù)，受傷人數(shù)序列各階自相關(guān)系數(shù)為：r1=0.8633，r2=0.6697，r3=0.4693，r4=0.2874，r5=0.0665。自相關(guān)系數(shù)歸一平均化，作為各滯時(shí)馬爾可夫鏈權(quán)重：w1=0.3664，w2=0.2842，w3=0.1992，w4=0.1220，w5=0.0282。

2.4 不同步長(zhǎng)轉(zhuǎn)移概率矩陣

根據(jù)2010年1季度至2019年第二季度共38個(gè)季度工傷數(shù)據(jù)，馬爾可夫工傷序列統(tǒng)計(jì)得不同滯時(shí)轉(zhuǎn)移概率矩陣：

2.5 創(chuàng)建預(yù)測(cè)表

以2018年第二季度至2019年第二季度共5個(gè)季度工傷人數(shù)、狀態(tài)轉(zhuǎn)移概率矩陣對(duì)2019年第三季度工傷人數(shù)進(jìn)行預(yù)測(cè)分析（見(jiàn)表4）。

表4 2019年第3季度工傷事故人數(shù)預(yù)測(cè)表

2019 年第三季度實(shí)際工傷人數(shù)為24 人/季度，由表可知，max(pi,i∈E)=0.5776，i=3，即預(yù)測(cè)2021年第三季度工傷人數(shù)狀態(tài)為3，工傷人數(shù)x滿足：26≤x＜39，預(yù)測(cè)數(shù)值與實(shí)際受傷人數(shù)相近，預(yù)測(cè)有效。取調(diào)整因子β=1，得狀態(tài)特征值μ=3.131。下限法得2019年第三季度工傷人數(shù)預(yù)測(cè)值：

x=+此值與實(shí)際相對(duì)誤差僅有1.25%，以2010年第1季度至2019年第3季度共39個(gè)季度工傷人數(shù)序列預(yù)測(cè)2019年第4季度的工傷狀態(tài)，得各階自相關(guān)系數(shù)、各步長(zhǎng)馬爾可夫鏈權(quán)重（見(jiàn)表5）。

表5 2010.1-2019.3季度序列各階自相關(guān)系數(shù)及權(quán)重

2010年第一季度至2019年第3季度共39個(gè)季度工傷事故受傷人數(shù)，不同步長(zhǎng)（滯時(shí)）轉(zhuǎn)移概率矩陣：

以2018年3季度至2019年第3季度共5個(gè)季度工傷人數(shù)預(yù)測(cè)2019年第4季度的工傷人數(shù)，預(yù)測(cè)見(jiàn)表6。

表6 2019.4季度工傷事故人數(shù)預(yù)測(cè)

2019 年 第4 季 度 實(shí) 際 工 傷 人 數(shù) 為22 人/ 季 度，max(pi,i∈E)=0.4248，i=2，即預(yù)測(cè)2019年第4季度工傷人數(shù)狀態(tài)為2，受傷人數(shù)x滿足：13≤x＜26，預(yù)測(cè)結(jié)果與實(shí)際情況完全一致，預(yù)測(cè)有效。

取調(diào)整因子β=0.95，狀態(tài)特征值μ=2.0231。下限法得2019年第4季度人工傷預(yù)測(cè)值：

此值與實(shí)際相對(duì)誤差僅有1.36%

2.6 不同步長(zhǎng)馬爾可夫鏈特征分析

此鏈狀態(tài)關(guān)聯(lián)互通，對(duì)任意i,j∈E,i←→j,(i≠j)，非周期，故為不可約正常返鏈。設(shè)平穩(wěn)分布為(πj,j∈E)：以步長(zhǎng)為2的轉(zhuǎn)移概率矩陣p(2)，平穩(wěn)分布、極限分布、各狀態(tài)重現(xiàn)期見(jiàn)表7。

表7 平穩(wěn)分布、極限分布及各狀態(tài)復(fù)現(xiàn)期

狀態(tài)j的復(fù)現(xiàn)期為T(mén)j，其對(duì)應(yīng)概率pj=1/Tj=πj。各狀態(tài)的復(fù)現(xiàn)期為T(mén)1=19.2678（季度），T2=3.8536（季度），T3=2.7525（季度），T4=5.6721（季度），T5=6.2696（季度）。在中長(zhǎng)期預(yù)測(cè)中，企業(yè)工傷人數(shù)為狀態(tài)3的概率最大，平均每2.7525季度重現(xiàn)，復(fù)現(xiàn)概率為0.3633；其次復(fù)現(xiàn)的為狀態(tài)2，每3.8536季度復(fù)現(xiàn)1次，概率為0.2595；復(fù)現(xiàn)概率最低的為狀態(tài)1為0.0519，復(fù)現(xiàn)期19.2678季度。

3 結(jié)論

1）加權(quán)馬爾可夫模型狀態(tài)的復(fù)現(xiàn)期T1=19.2678（季度），T2=3.8536（季度），T3=2.7525（季度），T4=5.6721（季度），T5=6.2696（季度）。長(zhǎng)期預(yù)測(cè)中，企業(yè)工傷人數(shù)狀態(tài)3的概率最大，平均每2.7525季度重現(xiàn)，復(fù)現(xiàn)概率0.3633；其次復(fù)現(xiàn)的為狀態(tài)2，每3.8536季度復(fù)現(xiàn)1次，概率為0.2595；復(fù)現(xiàn)概率最低的為狀態(tài)1，為0.0519，復(fù)現(xiàn)期19.2678季度。

2）實(shí)際應(yīng)用價(jià)值在于以某城市溫室氣體排放序列來(lái)預(yù)測(cè)未來(lái)溫室氣體年排放狀況，以某國(guó)家煤炭資源消耗量序列來(lái)預(yù)測(cè)未來(lái)煤炭年消耗狀況等。