---

2019 年 5 月

第 34 卷，第 5 期

[機器學習服務]

使用生存分析進行預測性維護

藉由Zvi Topol |2019 年

幾年前，我介紹存活分析的基本概念，並說明如何實作非參數化的演算法中呼叫 Kaplan Meier C# (msdn.com/magazine/dn630650)。現在，我打算探討另一部存活分析，特別是在兩個更進階的方法，都能在兩個常用的機器學習服務平台，Spark 機器學習程式庫 (MLLib) 和h2o.ai，兩者所支援的 Azure HDInsight。我將使用的預測性維護使用案例，在進行中範例所示。

工業物聯網的預測性維護

背後工業物聯網 (IIoT) 的主要概念是連接的電腦、 裝置、 感應器和工業設備與組織內的應用程式，並以持續收集資料，例如系統錯誤和電腦的遙測，從所有這些項目與分析，並根據這項資料，以獲得最佳的運作效率的目的。

預測性維護的目標是要準確地預測時的機器，或其中任何元件將會失敗。如果您可以這樣做，您可以執行維護，預測這類失敗發生之前。這是比未執行任何維護，直到發生失敗，在此情況下，電腦或元件，才可以使用失敗固定的如果確實是可修復更有效率。這類未規劃的停機時間很可能會很高。

預測性維護也會比在頻繁的間隔，也可能是來因為不必要的維護可能套用執行預防性維護更有效率的。

範例，我將使用的資料會在此範例是適用於的舊版bit.ly/2J4WnbN。此範例也包含製造 100 部的機器，與機器之間的任何相互依存性。每一部機器可以是其中一個可能的四個模型。

機器的資料會包含失敗、 維護作業和感應器的遙測，以及模型和存在時間 （以年） 機器的相關資訊的歷程記錄。此資料可在.csv 檔案可從先前所述的資源。我也會提供已轉換的資料檔案 (comp1_df.csv) 是 「 求生分析已準備好 」，並將說明如何在稍後執行的轉換。

原始範例中的每部電腦有四個不同的元件，但是我要把重點放只能在一個元件上。元件可以主動維護之前失敗，或進行修復的失敗後的維護。

存活分析

在存活分析有關我前一篇文章，我介紹了重要的基本概念，我要使用並擴充這篇文章中。建議您閱讀該文章，以熟悉這些概念，包括未回收和危險的函式、 censoring 和非參數化的 Kaplan Meier （公里） 估計工具。

在本文中，我將示範如何擴充公里估計工具，以包含 covariates 或變數 （也稱為功能） 可能會造成在存活，或是，在此情況下，機器元件的失敗影響的概念。在範例中，我將為 covariates 使用機器模型、 機器年齡和機器遙測，並使用求生迴歸模型來評估這類 covariates 對機器失敗。 

估計 covariates 對目標變數的概念在此情況下失敗、 危險速率或求生機率的時間、 不是唯一的存活分析和一般是迴歸模型的基礎。

當建立統計模型，您會看到 covariates 的三種主要的資料類型： 類別、 序數和連續。類別的資料類型所指分成幾個離散類別的類型。在這裡，機器模型是類別的資料類型，有四個不同的機器模型。序數的資料類型是類別的資料類型具有一些有意義的順序。比方說，從為 10，10 是最具娛樂性的其中一個，另一個的影片的評等最少。最後，連續資料類型是表示連續數字。這些就是機器遙測數據，也就是連續的數字 （在此情況下，每小時），在特定時間取樣。

之後用來識別要使用的資料類型和方法，您應該將各種資料類型編碼為 covariates。一般而言，迴歸模型的連續變數會自然地編碼為連續 covariates，而類別的資料類型將會需要某種形式的編碼方式。熱門選項這類編碼，我將在本文中，這是其中，含 N 類別分類的資料類型，會建立 N-1 covariates，和 i 的類別由將設定其特定共值之一，所有其他項目為零。第 n 個類別被以將所有 covariates 都設定為零。這通常是適合的迴歸模型，使用明確定義的基準，其中所有 covariates 可以都是等於零。這也是 R 程式設計語言使用來編碼類別變數或因素的格式。

Categoricals 這個編碼方式已設為零的部分或所有 covariates 的意義直接解譯。不過，連續資料類型，將特定共設定為零不一定有意義。例如，如果共代表機器高度或寬度，設定為零，共會是沒有意義，因為實際上沒有任何這類機器。

一種方式解決這個問題是使用平均值置中的連續 covariates，其中針對指定的共，其平均值表示定型資料上的會減去該值。然後，當您設定該轉換的共為零，它就相當於原始共設為其平均值。這項技術稱為 「 表示置 」，我將它用於機器年齡和遙測 covariates。

請務必記住，遵循此轉換，您應該一律使用平均值置中對齊的 covariates 做為模型的輸入。這也是如此時將套用至新的測試資料集的迴歸模型。

之後的資料值會編碼為 covariates，求生迴歸模型然後採用這些 covariates 和某些形式的存活目標變數 （我將討論推出），並指定模型繫結在求生/時間-至-事件這類 covariates 的影響.

轉換至求生格式和特徵工程設計資料

若要使用我將說明求生迴歸模型，您的資料必須要有至少兩個欄位： 感興趣 （這裡，機器故障） 和布林值欄位，指出是否 censoring 發生事件的時間戳記。（在這裡，censoring 描述或指定的時間之前，發生任何失敗的情況。在我的範例中，預防性的方式，而不做為回應失敗，發生的維護會視為 censoring。

我將討論求生迴歸模型具有不同來簡化其數學衍生所做的假設。部分這些假設可能不在這裡，保存，但仍可套用此範例模型的存活。

求生分析文獻是非常豐富許多進階求生迴歸模型及技術開發的位址和放寬這些假設部份。您可以深入了解這類模型中與技術的書籍，< 統計分析失敗時間資料 > Kalbfleisch 和 Prentice (Wiley-Interscience，2002年)，在bit.ly/2TACdLR。

我會假設每個機器元件上完全執行的維護作業會重設該元件，因此可以單獨處理。就可以使用兩種類型的間隔的存活迴歸 (示**[圖 1**):

• 失敗和上述的維護作業 （事件時間） 之間的間隔。
• 後續的維護作業 (censoring) 之間的間隔。

[圖 1 求生的表示法的機器失敗

在每個間隔**[圖 1**開始維護作業。間隔結束 X，其表示失敗，而第二種結尾 O，表示另一個維護作業之前失敗 （這是本質上的主動維護作業），這在此情況下表示 censored 的觀察第一個型別。

因此，原始的資料必須轉換成此格式包含兩個必要欄位。「 Time_to_event"欄位代表小時，直到發生失敗或下一個維護時間。「 事件 」 欄位設定為其中一個失敗，並維護作業失敗前為零。

最好經常在 covariates，執行額外的轉換通常稱為 「 特徵設計 」。 此程序的目的是產生 covariates 與更好的預測功效。例如，您可以建立另一個會計算平均值的壓力失敗前 10 個小時內的共。有許多不同的選項函式和可能的時間範圍來建立這類 covariates，還有一些工具可用來協助自動化此程序，例如開放原始碼 Python 封裝 tsfresh (tsfresh.readthedocs.io/en/最新)。

現在我要討論的兩個求生迴歸模型： Cox 成比例的危險模型 （或 Cox PH 模型） 提供 h2o.ai 和 Spark MLLib 中的 Weibull 加速失敗時間模型。

Cox 成比例的危險迴歸

您應該記得危險函式，會決定事件速率，在時間 t 物件或個人存留在時間 t。預測性維護的範例中，它可以描述為下一個小時的時間，針對給定的時間 t 和自其上次維護元件 1 失敗尚未發生的所有機器失敗的機率。更高的危險率表示較高的風險發生失敗。Cox PH 迴歸評估下列模型所指定的危險速率 covariates 的影響：

在這裡，h(t) 是危險函式，在時間 t h0(t) 基準危險在時間 t、 Xi 變數都是不同的 covariates，是對應的 beta 版會對應至 covariates 的係數 (有更多一點的更新版本)。基準危險時所有 covariates 都會等於零的危險。請注意，這與密切相關的截距，在其他的迴歸模型，例如線性或羅吉斯迴歸。

根據此模型中，covariates 的存活時間之間沒有直接關聯性。此模型稱為半參數化，因為在時間 t 的危險速率是基準危險的一項速率，估計的資料，並不需要參數化的已關閉的表單並參數化的乘法類元件的函式。

此模型稱為調和間距的危險模型的原因是因為它可讓您比較兩個障礙函式的比率。在此情況下，提供預估的模型，是兩個不同的資料點之間的比率：

基準危險速率抵銷功能和所處理的風險之間產生的比例是係數和 covariates 函式，並且一次未相依於時間。這是密切相關羅吉斯迴歸估計機率的記錄檔的位置。此外，Cox PH 迴歸模型並不直接指定求生函式，並同時提供焦點比率或危險函式的比例的資訊。因此，它主要用來了解 covariates 對生存能力，而不是直接估計求生函式。

Cox PH 與指定的模型、 係數和非參數化的基準危險可估計使用各種技術。一個常用的技巧是部分的最大可能性估計 （也用於 h2o.ai）。

下列程式碼片段是執行估計的平均值置中對齊 covariates （機器遙測和年齡） 上使用 h2o.ai Cox PH 模型類別共機器模型的 R 指令碼：

library(h2o)
localH2O <- h2o.init()
inputFileName<-'comp1_df.csv'
df<-read.csv(inputFileName, header=TRUE, stringsAsFactors=TRUE)
df.hex <- as.h2o(df, key = "df.hex")
model <- h2o.coxph(x = c("age_mean_centered", "model","volt_mean_centered",
  "rotate_mean_centered","pressure_mean_centered", "vibration_mean_centered"),
  event_column = "event", stop_column ="time_to_event" ,training_frame = df.hex)
summary(model)

在撰寫本文時，h2o.ai Cox PH 模型無法使用從 Python，因此會提供 R 程式碼。安裝指示位於bit.ly/2z2QweL，或針對使用 Azure HDInsight 的 h2o.ai 在bit.ly/2J7nXp6。

執行程式碼片段會產生輸出所示**[圖 2**。

[圖 2 輸出 Cox PH 迴歸

Surv(time_to_event, event) ~ model + volt_mean_centered + rotate_mean_centered +
    pressure_mean_centered + vibration_mean_centered + age_mean_centered
  n= 709, number of events= 192
                               coef      exp(coef) se(coef)    z      Pr(>|z|)   
model.model2                  -0.066955  0.935237  0.257424 -0.260    0.795   
model.model3                  -0.021837  0.978400  0.215614 -0.101    0.919   
model.model4                   0.308878  1.361896  0.227469  1.358    0.174   
volt_mean_centered             0.031903  1.032418  0.003990  7.995    1.33e-15 ***
rotate_mean_centered           0.001632  1.001633  0.001362  1.199    0.231   
pressure_mean_centered        -0.008164  0.991869  0.005768 -1.415    0.157   
vibration_mean_centered        0.018220  1.018387  0.013866  1.314    0.189   
age_mean_centered              0.004804  1.004815  0.013293  0.361    0.718   
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
                                  exp(coef) exp(-coef) lower .95 upper .95
model.model2                      0.9352     1.0692    0.5647     1.549
model.model3                      0.9784     1.0221    0.6412     1.493
model.model4                      1.3619     0.7343    0.8720     2.127
volt_mean_centered                1.0324     0.9686    1.0244     1.041
rotate_mean_centered              1.0016     0.9984    0.9990     1.004
pressure_mean_centered            0.9919     1.0082    0.9807     1.003
vibration_mean_centered           1.0184     0.9819    0.9911     1.046
age_mean_centered                 1.0048     0.9952    0.9790     1.031
Rsquare= 0.094   (max possible= 0.941 )
Likelihood ratio test= 70.1  on 8 df,   p=4.69e-12
Wald test            = 70.19  on 8 df,   p=4.514e-12

請注意第一個重點是 covariates 估計的係數。機器模型共會編碼成類別資料型別。此類別的基準是 model1 表示 splittunneling 編碼方式的其他三個機器模型 （model.model2、 model.model3 和 model.model4） 為零的三個 covariates。每個共取得自己的係數。了解如何解譯係數相當重要的。

如果您對機器模型 covariates (exp(coeff) 在輸出中)，您會看到該 model.model2 有 model.model4 有值為 1.3619 0.9352，值的係數套用指數函式。這表示 model2 的機器具有危險速率低於基準機器模型 （模型 1）、 危險率的 6.5 %model.model4 機器有相較於機器 model.model1 36.2%的大幅上升危險。換句話說，model.model4 的機器都有最高失敗的風險，而機器 model.model2 具有最低的失敗風險。因此，當排列優先順序的維護作業，機器的模型應納入考量的重要因素。

所有其他 covariates 是 mean 置中的連續 covariates。與其建立關聯的係數的解譯是現在的危險比例由所提供的指數，其方法周圍 covariates 也。因此，增加一個單位 （保留所有其他固定的 covariates） 共值，危險比率增加 （或減少） 的指數的係數 (類似的方式的類別變數)。因此，比方說，藉由增加一個單位電壓，失敗的風險會增加百分之 3.2。

此處提及的另一個重點是有關模型診斷技術。這種技術提供了解是否適當 （在此情況下，Cox PH 模型） 被視為模型的基礎。在這裡，Rsquare 值 （介於零和一段，愈高愈好之間的值） 是相對較低 (0.094) 和大部分的係數 z 分數不表示係數是統計上明顯 （沒有足夠的辨識項，以便支援它們不同於零）。這兩個這些指標會導致沒有改進，例如透過特性工程設計空間的結論。另外還有其他專屬於 Cox PH 模型，應該進行的統計檢定。如需詳細資訊，可以稍早查閱求生分析文獻中所述。

加速的失敗時間模型

Spark MLLib 中的存活迴歸模型是加速失敗時間 (AFT) 模型。此模型中直接指定某些理論的數學發佈 (Weibull) 求生函式，並加速的失敗時間屬性。

SITE 的模型定義，如下所示。假設物件以使用 （線性） covariates 和係數：

也假設物件已參數化生存函式 s(t) 並，以 s0(t)，基準 （具有物件設定為零的所有 covariates) 求生函式。SITE 的模型會定義 s(t) 為 s0(t) 之間的關聯性：

從這個定義中，您可以看到模型為何稱為加速失敗時間模型。這是因為求生函式包含加速器因數，也就是線性組合的 covariates，乘以存活時間 t 指數函式。

某些 covariates，例如使用期限時，這種模型會很有用 （在 [我的資料集，電腦年齡），可能會造成單純加速或減速的存活/失敗的時間。

韋伯分佈是指數分佈的概括，常用參數化生存模型中的連續分佈。有一些有關如何將它參數化的變化形式。在這裡，我將使用下列兩個參數 Weibull 分佈版本的 t > = 0:

（另外還有三個參數的版本。） 分佈的兩個參數是取決於 k 和級別，取決於 lambda 的圖形。粗略的比喻是鐘型分佈的特性的平均值和標準差的方法。

回想一下，發佈密度函式 f(t)、 危險函式 h(t) 和求生函式 s(t) 之間的關聯性根據 f(t) = h(t)s(t)。

韋伯危險和求生函式如下：

不同 Cox PH 模型中，於生存和危險函式完整指定並有參數化的表示法。請參閱**[圖 3並[圖 4**的 k 和 lambda 的不同值的 Weibull 分佈和求生函式的視覺效果。

[圖 3 Weibull 分佈圖形為不同值的 K 和 Lambda 的函式

[圖 4 Weibull 求生的 K 和 Lambda 的不同值的函式圖形

[圖 5說明 SITE 模型 covariates 有 Weibull 求生函式的圖形的影響。

[圖 5 加速失敗時間 Weibull 求生機率函式

Spark MLLib 中 AFT Weibull 模型的係數的估計是使用最大可能性估計演算法。您可以深入了解如何在完成bit.ly/2XSauom，並尋找實作程式碼，在bit.ly/2HtJw0v。

不同於這項估計 Cox PH 模型，其中的 covariates 係數回報 （以及一些診斷），結果會取自估計 Weibull AFT 模型報表 covariates，以及特定參數的係數韋伯分佈的 — 攔截和小數位數參數。我將示範如何將轉換的 k 和位元中的 lambda。

韋伯 AFT 實作 Spark MLLib 中的結果符合 Weibull AFT 實作使用熱門的 R 程式庫 」 求生"survreg 函式的結果 (更多詳細資料位於bit.ly/2XSxkw8)。

您可以執行下列 R 指令碼的 AFT Weibull 模型估計 (在本機安裝的 Spark MLLi 上, 執行的程式碼，但您也可以在 HDInsight 上使用 Spark bit.ly/2u2U5Qf):

library(survival)
library(SparkR, lib.loc = c(file.path(Sys.getenv("SPARK_HOME"), "R", "lib")))
sparkR.session(master = "local[*]")
inputFileName<-'comp1_df.csv'
df<-read.csv(inputFileName, header=TRUE, stringsAsFactors=TRUE)
aftMachineDF<-suppressWarnings(createDataFrame(df))
aftMachineModel <- spark.survreg(aftMachineDF,Surv(time_to_event, event) ~ model +
  age_mean_centered +  volt_mean_centered + rotate_mean_centered +
  pressure_mean_centered + vibration_mean_centered)
summary(aftMachineModel)

指令碼將產生的估計的係數含其他資訊。就能夠執行 survreg （因為符合的結果），以取得這類資訊：

library(survival)
machineModel<-survreg(Surv(time_to_event, event) ~ model + age_mean_centered +
  volt_mean_centered+rotate_mean_centered + pressure_mean_centered +
  vibration_mean_centered, df, dist='weibull')
summary(machineModel)

在此情況下，R 指令碼會產生更複雜的輸出所示**[圖 6**。

[圖 6 輸出 Weibull SITE 迴歸

survreg(formula = Surv(time_to_event, event) ~ model + age_mean_centered +
    volt_mean_centered + rotate_mean_centered + pressure_mean_centered +
    vibration_mean_centered, data = df, dist = "weibull")
                              Value        Std. Error    z            p
(Intercept)                 8.172991       0.119133    68.6040     0.00e+00
modelmodel2                 0.040289       0.154668    0.2605      7.94e-01
modelmodel3                 0.027225       0.129629    0.2100      8.34e-01
modelmodel4                  -0.163865      0.136382   -1.2015     2.30e-01
age_mean_centered           -0.000753      0.007960    -0.0946     9.25e-01
volt_mean_centered           -0.019731      0.002583   -7.6391     2.19e-14
rotate_mean_centered         -0.000767      0.000821   -0.9334     3.51e-01
pressure_mean_centered       0.005173       0.003496     .4795     1.39e-01
vibration_mean_centered      -0.008214      0.008391   -0.9789     3.28e-01
Log(scale)                   -0.508060      0.051963   -9.7773     1.41e-22
Scale= 0.602
Weibull distribution
Loglik(model)= -1710.3   Loglik(intercept only)= -1747.2
                           Chisq= 73.73 on 8 degrees of freedom, p= 8.9e-13
Number of Newton-Raphson Iterations: 8
n= 709

在我們繼續之前描述的輸出中，我應該提醒 Weibull 參數化在 Spark MLLib 和 survreg 是有點不同於我所討論的參數化。

轉換需要，而且可以進行，如下所示。代表參數報告 — 截距的 m 和小數位數 s，然後 k = 1/秒，lambda = exp(-m/s) 和每個應該乘以係數 (-1/s)。沒有 R 封裝呼叫可以自動執行這項轉換的 SurvRegCensCov 模型上使用 ConvertWeibull 該 survreg 估計：

$vars
                                            Estimate        SE
lambda                                    1.260459e-06   8.642772e-07
gamma                                     1.662064e+00   8.636644e-02
modelmodel2                               -6.696297e-02  2.569595e-01
modelmodel3                               -4.524990e-02  2.155000e-01
modelmodel4                               2.723541e-01   2.268785e-01
age_mean_centered                         1.251958e-03   1.322780e-02
volt_mean_centered                        3.279500e-02   3.947495e-03
rotate_mean_centered                      1.274045e-03   1.365339e-03
pressure_mean_centered                    -8.598142e-03  5.807130e-03
vibration_mean_centered                   .365213e-02    1.391255e-02

在這裡，gamma 等於 k 從先前的 Weibull 參數化。(如需有關 SurvRegCensCov 的詳細資訊，請參閱bit.ly/2CgcSMg。)

估計的參數，指定與不同的是，Cox PH 模型中，現在正是可以直接取得求生函式 （它是 Weibull AFT 求生函式），並使用它來預測任何 covariates 求生機率。假設資料集內的第一個點新的資料點，您可以執行下列命令：

predict(machineModel, newdata=df[1,], type='quantile')

這會產生的時間 （以小時為單位） 的事件分位數 0.1 和 0.9 （預設值），就像這樣：

807.967 5168.231

這表示，假設 （此處列出） 的第一個資料點的 covariates，故障的機率是百分之 10 在或之前 807.967 維護作業，之後的小時，而故障的機率是在或之前 5168.231 小時的 90%下列維護作業：

model       age           volt_mean_centered                 rotate_mean_centered
 model3       18                  3.322762                          51.8113
pressure_mean_centered     vibration_mean_centered             age_mean_centered
 10.10773                         11.4267                           6.488011

您也可以使用參數"p"以取得零和一; 之間的任何分位數的存活時間例如，將參數新增 「 p = 0.5 」 可讓中間失敗時，第一個資料點，這是 2509.814 維護作業之後的小時。

如同 Cox PH 模型評估中，p 中的資料行的 survreg 輸出會提供預估，不過在此情況下將數字的較佳的係數 （較低 p-值） 的統計重要性相關資訊。仍有功能工程這裡所 Cox PH 模型之前所述的空間。

務必也執行模型的診斷，如同 Cox PH 迴歸，並確定 Weibull AFT 模型是適合的資料，相較，比方說，其他的參數化模型中。雖然我將不會說明此程序，您可以深入了解所指的 [求生分析] 活頁簿先前所述。

總結

我在提供用於預測性維護 IIoT 為採用兩個求生迴歸模型中所提供的激勵範例h2o.ai和 Spark MLLib。我示範了如何建立求生分析 framework 中的機器失敗預測性維護問題模型 covariates 和轉換時間序列資料為求生格式編碼的變數。我也說明了求生兩個模型，以及如何將它們套用至資料之間的差異。最後，我談到簡短的結果和模型診斷的解譯。請務必請注意，我只不過粗略介紹本有趣且非常豐富的主題，而且我建議您深入探索。藉由參考我的文章中提及的文宣因此是這麼做的起始點。

---

Zvi Topol開始從事資料科學家在各種產業有關，包括行銷分析、 媒體和娛樂和工業物聯網。他已傳遞，並會導致多個機器學習服務與分析專案，包含自然語言、 語音介面、 認知服務，影片分析、 推薦系統和行銷決策支援系統。Topol 正 MuyVentive llc 的進階的分析 R & D 公司，可以在達到zvi.topol@muyventive.com。

感謝下列 Microsoft 技術專家來檢閱這篇文章：James McCaffrey

---

MSDN Magazine 論壇中的這篇文章的討論