基于Rényi熵的q-指數分布及其可靠性分析應用

謝 暄，王敏夷，白穎利，汪東敏，李西峰*，謝永樂

(1. 電子科技大學自動化工程學院 成都 611731；2. 中國空間技術研究院通信衛星事業部 北京 海淀區 100094；3. 四川慧龍科技有限責任公司 成都 610041)

隨著集成系統復雜性的日益增加，系統建模和可靠性保證面臨越來越大的挑戰。在可靠性分析的各種場景中，雖然指數分布模型已被較廣泛地用于建立系統壽命模型[1]，但由于指數分布的風險函數是恒定的，因此直接將指數分布用于描述系統壽命，存在無法描述損傷過程和無法準確反映故障累積效果的問題，最典型的例子是將指數分布簡單應用于描述人類死亡率和電子設備生命周期，效果不夠理想。對這類過程，通常需要采用具備浴盆特征的風險函數所對應的壽命分布來準確描述。韋伯分布作為指數分布的概括，以及它的帶有浴盆型風險函數的擴展得到了重視，并被廣泛應用于許多領域[2-4]。

目前關于指數分布推廣的研究，主要集中在廣義韋伯分布的累積分布函數的參數加法上。大多數情況下，這些韋伯類型的一般化是通過簡單的參數加法技術來實現的。此方法一般除保留了原有的參數，還引入了一些新的參數，使得模型對壽命數據的擬合，多數情況下優于沒有新參數的模型。

基于參數加法，文獻[5]提出了四參數廣義指數模型，稱為廣義指數化線性指數分布。結果表明，這一擴展可以導出一系列指數型分布，如指數化韋伯分布。結合韋伯分布和改進的五參數韋伯分布，文獻[6]提出了一種修正的韋伯分布擴展?；诖?，文獻[7]又提出了一個離散五參數修正型的韋伯分布，并發現基于這種分布的離散數據擬合勝過其他3個修正的韋伯模型。

盡管韋伯推廣在壽命模型中具有良好的效果，但由于參數估計過程復雜，參數物理意義不明確(至少這些參數在這些廣義的分布中，不能直接指示系統的當前狀態)，限制了韋伯推廣的應用[8-10]。

本文提出了一種廣義雙參數指數分布，它可以作為可靠性分析的基礎，例如浴盆型風險函數的構造和有用壽命預測。與傳統的參數加法不同，本文用最大熵原理得到廣義指數分布。從物理角度看，這種廣義分布的主要優點是兩個參數都具有明顯的物理意義。一個參數具有分形意義，表示系統的穩定程度；另一個參數則與系統的平均行為密切相關。兩個參數均可以用來刻畫系統的性能。

1 q-指數分布模型

1.1 一般化的指數分布

基于Rényi熵，可以通過最大熵方法，得到廣義指數分布：

式中，q是分形參數；λ是分布的期望值；q和λ都是非負的。根據參數q和λ的值，概率密度函數fq會有不同的形態。圖1所示為不同的參數q時固定均值λ下的fq形態?？梢钥闯鰍-指數分布是單峰分布，也是右偏分布。因此，q-指數分布可以被考慮用于對壽命相關的可靠性分析。

圖1 q-指數分布fq(x)

計算可得，q-指數分布的累積分布函數是：

式中，2F1是高斯超幾何函數。

1.2 特殊衍生分布

下面將給出q-指數分布的兩個主要的衍生分布，它們在生存分析、壓縮感知、剩余使用壽命預測和其他的可靠性相關領域有潛在實用價值。

首先，當q→1時，該分布簡化為均值為λ的指數分布：

其次，使用式(1)可以得到一個廣義的q-拉普拉斯分布：

圖2顯示了不同參數值q和固定λ值(λ=10)的廣義拉普拉斯分布。拉普拉斯分布在壓縮感知領域起到了重要作用[11-13]，這里衍生出的拉普拉斯分布在統計信號處理和機器學習領域具有潛在應用價值。

圖2 不同參數q下的q-拉普拉斯分布gq(x)

2 模型性質

2.1 q-指數分布的風險函數

生存函數是表示一系列事件的隨機變量函數，通常用于表示一些基于時間的系統失敗或死亡概率。假設T表示產品使用壽命，其分布函數為F(t)，那么該產品壽命大于t的概率為：

S(t)=P(T＞t)=1-F(t)

式中，S(t)被稱為生存函數。在此基礎上，可用風險函數刻畫已有效使用到t時刻的產品，在[t,t+Δt]極短時間內“死亡”的風險：

據此計算可得q-指數分布的生存函數：

風險函數：

圖3為具有不同參數值的風險函數圖像，可以看到風險函數具有許多不同的形狀，q-指數分布中當q＞1，顯示遞增特性；當q＜1，顯示遞減性。

圖3 q-指數分布的風險函數

注意到q-指數分布的風險函數呈現多態性，其中部分具有澡盆特征?？筛鶕嶋H使用環境，選擇恰當的q-指數建模風險過程，從而可以提高壽命估計精度。

2.2 q-指數分布的矩

q-指數分布fq的均值是λ，方差是：

可以得到分布的k階矩為：

式中，k表示分布的k階矩。

3 參數估計

為了在不同場景下選擇恰當的參數，這里給出兩種估計q-指數分布中參數的方法：最大似然估計法(maximum likelihood estimation,MLE)和信息似然估計法。最大似然估計法適用于滿足高斯分布的數據集，信息似然估計法適用于高斯和非高斯分布的數據集。

3.1 最大似然估計

根據式(1)，建立對數似然函數：

然后有：

q和λ可以通過求解方程組(6)獲得：

3.2 信息似然估計

另一個估計Rényi信息未知參數的方法，將信息論和譜估計相結合，可在最小先驗的條件下，取得經驗風險最小的參數估計值[14]。具體步驟如下：

首先，已知Rényi信息頻譜定義為：

當λ →1，

Rényi信息頻譜趨向于香農熵。

同時，Rényi信息的頻譜梯度被定義為：

式中，D(·)是隨機變量x的方差。

令φf:=-2LR(1)

一方面，將式(1)中的fq代 入式(9)計算φf，根據文獻[13]可得:

另一方面，使用核方法估計 φf。具體而言，φf可以通過核方法計算如下：

式中，

fn(x)是為無參數內核密度估計量，定義為：

假設K(·)為有界變分的概率密度函數(核函數)，其支撐集位于部分有限區間。Xi是隨機變量X的樣本。設bn為滿足以下收斂條件的序列：

比較式(10)和式(11)，式(1)中的參數q可以通過求解式(13)得到：

4 數據分析

本文使用了3個與可靠性相關的實驗去評估q-指數分布在可靠性分析中的有效性，包括白血病病人生存期分析、設備壽命預測以及鋰電池壽命預測。

第一個數據集是40個白血病病人的生存期數據[15]，它具有不斷增加的風險率。第二個數據集來自50個組件的樣本[16]，具有浴盆型風險的風險率。利用K-S統計量評估了q-指數與上述兩個著名數據集的擬合結果，并計算出赤池信息準則(AIC)值和貝葉斯信息準則(BIC)值，用于比較它和其他廣義指數分布的擬合優越度。第三個實驗建立了關于3組鋰電池的容量退化模型，與普通指數分布的結果相比，利用q-指數分布可以得到更準確的剩余壽命預測結果。

4.1 白血病病人生存期

該數據集由文獻[15]給出，數據源自沙特阿拉伯衛生醫院部門，如表1所示。它記錄了40名白血病患者的生命周期。

表1 40名白血病患者的生存期

基于雙參數模型的生存函數經驗估計如圖4所示，圖5描述了模型的經驗和擬合風險函數。同時，由于該數據集所示患者死亡風險是隨時間上升的，因此，韋伯分布、線性失效率分布(linear failure rate distribution, LFR)和q-指數分布均可作為數據擬合的候選者。為了判斷上述哪個密度函數更適合數據擬合，用MLE方法對3種模型參數進行了估計，如表2所示。然后利用幾種不同的測試統計測度對擬合結果進行了評估。

表2 白血病數據集各模型參數的MLE值

圖4 白血病數據集的生存函數

由統計學可知，K-S值表示分隔程度，一般大于0.2即表明模型具有良好的分隔性能。由表3可知，q-指數分布與韋伯分布、指數分布和LFR分布的K-S均有良好分隔能力。同時，由表3列出的對數似然函數值可知，q-指數分布與韋伯分布、指數分布的最大似然函數值均在相近水平，所以q-指數分布可以充分利用先驗信息獲得對未知參數的最大似然估計。

表3 白血病數據的對數似然函數值和K-S

根據統計學原理，p值是判斷原假設是否成立的依據，一般認為p＞0.05，說明兩組樣本無統計學差異。由表4可知，q-指數分布與韋伯分布的p值在相近水平，q-指數分布具備描述壽命分布的能力。同時，通過計算AIC或者BIC值，相比指數分布、韋伯分布和LFR分布，不難發現q-指數分布具有更小的AIC或BIC值，因此q-指數分布具有更好的壽命數據擬合性。

表4 白血病數據集各模型的p值、AIC和BIC

4.2 元件壽命數據

如表5所示，這是由文獻[17]提供的50臺設備的壽命數據，已有研究人員利用韋伯分布[18-21]、指數分布[5]、LFR分布[22]分析了這個數據集。表6給出了所使用的每個分布參數的MLE估計值。圖6和圖7分別給出了設備數據集的經驗參數生存函數以及風險函數。

圖6 Aarset數據集的生存函數

圖7 Aarset數據集的風險函數

表5 50臺設備的生存期

表6 設備時間數據集各模型參數的MLE值

已知Aarset數據集的風險函數具有浴盆形狀。不失一般性和為了分析簡便，本文只專注風險函數中單調遞增部分，使用雙參數分布來近似風險函數。為了進行參數比較，使用似然檢驗去對比原假設和備選假設。此外，利用AIC[23]在多個模型中選擇最優模型。最適合數據擬合的模型應具有最低的AIC。表7和表8給出了對數似然函數值、KS值、p值、AIC和BIC值[24]。

表7 設備時間數據集各模型的對數似然函數值和K-S值

表8 設備時間數據的P值，AIC和BIC

對于Aarset數據集，由表7可知q-指數分布的K-S值大于0.2，所以具備良好的分隔能力。由所列對數似然函數值可知，q-指數分布與韋伯分布、指數分布的對數似然函數值在相近水平，說明q-指數分布亦可以充分利用先驗信息獲得對未知參數的最大似然估計。

由表8知，q-指數分布的p值大于0.05，因而具備描述Aarset數據集的能力。同時，可以看出q-指數分布在本文所提到的所有分布中具有最小的AIC和最小的BIC值。這說明在所列分布中，q-指數分布能夠最好地擬合本數據集。

4.3 鋰電池壽命預測

估計剩余使用壽命(remaining useful life,RUL)有助于降低實際系統中發生災難性事件的概率[25]。為了研究所提出的q-指數分布的有效性，本文采用美國宇航局(NASA) Prognostics公司(PCoE)的3個電池數據集來預測鋰電池的剩余使用壽命。在該數據集中用電池容量來表征壽命，容量越低，剩余壽命越少。

3組電池(即蓄電池005、蓄電池006和蓄電池018)屬于同一類型，通過在室溫下工作在3種不同的狀態下(充電、放電和阻抗)進行加速老化試驗[26]。這種電池的額定容量是2 Ah，當電池容量減少到額定容量的70%(從2 Ah減少到1.4 Ah)時，電池就會達到使用壽命終止(EOL)標準，容量數據如圖8所示。

圖8 NASA PCoE的電池容量數據

容量退化過程可用狀態空間模型來描述：

式中，xk表示k周期的真實容量值；yk表示k周期的預測值；wk-1表 示環境的擾動；υk表示觀測噪聲。

分別采用指數分布和q-指數分布作為h(xk)來描述容量狀態轉換，采用粒子濾波算法[27]自適應地確定預測容量值，并使用最大似然法確定指數分布和q-指數分布的相應參數。

由粒子濾波基本原理可知，按照頻率派的概率觀點[24]，這里剩余壽命預測的概率分布可以通過對粒子濾波每個周期間隔內的粒子個數計數得到。間隔內粒子數越多，說明該區間代表大多數情況下的壽命長度。

在電池005實驗中，共有168個循環樣本，分別用前60個、80個和100個樣本點訓練粒子濾波器，設定粒子濾波預測的起始時間分別為T=60、80、100。當起始點T=100時，RUL的分布直方圖和RUL預測結果如圖9、圖10和表9所示。其中，圖9b使用參數q=1.001。

圖9 電池005的RUL直方圖

此外，為了定量評價預測精度，將預測誤差定義如下：

式中，RULP表示預測的周期數；RULt表示真實的周期數。

如圖10(其中q=1.001)所示，電池005的真實壽命周期的結束數為124，使用q-指數分布相對能更好地擬合實際值。根據表9，在指數分布假設下，在起始預測點T=60時，粒子濾波算法預測的壽命周期數為152。根據上述定義，指數分布模型預測誤差為：

表9 不同起始點T下鋰電池的RUL預測結果

圖10 電池005在起始點T為100的RUL預測結果

eRUL=RULp-RULt=152-124=28

同樣以T=60為起始預測點，q-指數分布在q=1.001的情況下，預測誤差為23，小于指數分布假設的結果。因此，提出的q-指數分布輔助粒子濾波算法對電池005的剩余使用壽命有更準確的估計。

同理，觀察圖11(其中q=1.010)與圖12(其中q=0.990)中電池006和電池018的結果，結合表9，發現對于電池006而言，真實壽命周期的結束數為108，而在指數分布假設下，在起始預測點T=60、80、100時，使用q-指數分布獲得預測結果均與指數分布的相同，說明基于q-指數分布的粒子濾波方法至少可以獲得和指數分布假設下一樣的估計精度。

圖11 電池006在起始點T=100的RUL預測結果

圖12 電池018在起始點T = 80處的RUL預測結果

對于電池018而言，真實壽命周期的結束數為96，在指數分布假設下，在起始預測點T=60時，預測的壽命周期數為92，預測誤差為-4。而q-指數分布在q=0.990的情況下，如表9所示的同一點的預測誤差為-2，誤差預測減小了50%。對于T=80和T=100的起始預測點，使用q-指數分布均可取得更準確的效果。

5 結 束 語

本文基于最大熵方法，通過計算均值約束下最大Rényi熵，得到一種新的廣義指數分布：q-指數分布。本文對q-指數分布的統計特性進行了分析，并給出了均值和各階矩的解析公式。為了便于應用于可靠性分析，給出了基于q-指數分布可靠性的預測模型及對應的生存函數和風險函數的解析表達式。指出可采用了兩種方法：極大似然估計法和信息似然估計法進行雙參數估計。最后，結合醫學白血病患者壽命數據集、設備元件壽命數據集及鋰電池剩余壽命數據集進行了驗證，通過與韋伯分布、指數分布等常用壽命預測分布對比，驗證了q-指數分布的有效性和估計精度的優良性。下一步將挖掘q-指數分布在復雜系統建模中的高效應用。