[譯]利用貝葉斯推理做硬體故障率的準實時預測

原文連結（英文）： https://databricks.com/blog/2019/02/14/near-real-time-hardware-failure-rate-estimation-with-bayesian-reasoning.html 導語：你可能已經不知不覺中在資料科學項中用上了貝葉斯相關技術！如果你還沒用上，這個技術可以增強你的資料分析能力。本文會接着 Data Science Central 中這篇貝葉斯推理簡介文章，展示這項技術在現實世界中的應用案例：通過傳感器收集的流式資料預測硬體故障率。如果你對此感興趣，請繼續閱讀！

使用者案例：硬碟故障

如果你可以檢測資料中心的硬碟驅動器，擷取包括驅動器型号和故障等資訊。你應該想知道哪些硬碟型号的故障率最高。

Backblaze公司收集了近幾年内資料中心硬碟運作記錄，并公布了

它們的故障資料

。這些資料按單個驅動器和日期進行細分，并具有來自SMART傳感器讀數的豐富資訊，但本文将僅僅使用驅動器的型号以及它在某一天是否發生故障了的資訊。以2018年第二季度和第三季度資料為例，下載下傳CSV檔案的壓縮包并解壓縮到一個共享位置。

最直接的做法是按型号計算故障率，并将其作為實際故障率的最佳估計值。按故障率計算的前六類驅動器是（注意：此處的比例按年度計算為年的故障率，為便于閱讀）：

from pyspark.sql.functions import *

data_path = ...
agg_df = spark.read.\
  option("header", True).\
  option("inferSchema", True).\
  csv(data_path + "/current/").\
  select("model", "failure").\
  groupBy("model").\
  agg({"*": "count", "failure": "avg"}).\
  filter("avg(failure) > 0").\
  orderBy(desc("avg(failure)")).\
  limit(6).\
  cache()
display(agg_df.withColumn("pct_failure", col("avg(failure)") * 100.0))           

（在Databricks平台，隻需使用單元格輸出中的繪圖選項就可以将結果的表格視圖切換到上面的直方圖。） 型号 ST8000DM004 全年故障率是0.5%，看起來是一個很高的數值，但是，我們也無法确定它是否隻是在第二季度和第三季度經曆了高故障率而已。

假設有一個真實的故障率 λ，那麼機關時間内的故障數目應遵循泊松分布Poisson(λ)（就像抛硬币實驗中正面朝上的次數遵循二項分布）。 如上所述，我們希望能了解故障率λ本身的機率分布，它遵循一個帶兩個參數ɑ和β的伽馬分布（gamma distribution）。和貝塔分布（Beta distribution）一樣，他們有一個簡答的解釋：在t天内觀察到f次故障後，故障率λ的分布是Gamma(f,t)。

import numpy as np
import scipy.stats as stats
from pyspark.sql.types import *

to_gamma_schema = ArrayType(StructType([StructField("rates", DoubleType()), StructField("pdfs", DoubleType())]))

def build_gamma_udf(df, max_scale_factor=1):
  max_avg_failure = df.select(max("avg(failure)")).collect()[0][0]
  # Computes the PDF of the gamma distribution at points for each model's rate param
  def to_gamma_pdf(count, avgf):
    # k = alpha = number of failures
    k = count * avgf
    # theta = 1/beta, where beta = total observations
    theta = 1.0 / count
    # Pick rates to try with more samples at small rates
    rates = np.geomspace(max_avg_failure / 1000, max_avg_failure * max_scale_factor, 100)
    # Calculate the probability distribution function at those points
    pdfs = stats.gamma.pdf(rates, a=k, scale=theta)
    return list(zip(rates.tolist(), pdfs.tolist()))
  return udf(to_gamma_pdf, to_gamma_schema)

to_gamma_pdf_udf = build_gamma_udf(agg_df)

gamma_df = agg_df.\
  withColumn("rates_pdfs", explode(to_gamma_pdf_udf("count(1)", "avg(failure)"))).\
  select("model", col("rates_pdfs.rates").alias("rate"), col("rates_pdfs.pdfs").alias("pdf"))

display(gamma_df.withColumn("pct_rate", col("rate") * 100.0))           

我們注意到不同型号硬碟的預估故障率分布是非常不一樣的。型号ST8000DM004有着非常廣的機率分布，同時，不得不承認它也有着更高的故障率；為了顯示友善我們隻截取了部分曲線，但它的很大一部分機率分布在大于每年0.4%的區間。這是因為相對而言這個型号硬碟的觀察樣本比較少，見下表：

放大左側的兩個峰值分布，分别是“ST500LM012 HN”和“WDC WD5000LPVX”：

有趣的是，WDC可能性最高的年故障率是比ST500低的（大概是0.015% 比 0.02%）。但是，因為WDC的觀測樣本較少，是以它的可能故障率分布會更廣。如果我們的問題是哪種型号的故障率最有可能高于0.03%，基于目前的資訊，WDC反而是更有可能的。類似的場景比如說：我們需要決定哪個型号的故障最有可能超出制造商的可容忍上限。此時，不僅僅是最大後驗機率估計，機率分布也同樣重要。

用先驗機率分析資料流中的故障資訊

實際上，這些觀測資料并非一次性全部能擷取到。資料是以流的形式到達，是以很自然的我們需要持續的運作這些計算程式。 Apache Spark Structured Streaming非常适合用來處理這種資料，它可以幾乎完全相同的形式運作上述代碼。

當然，這類流式分析程式從第一天就開始運作時，并沒有任何曆史資料。甚至當積累了兩個季度的資料之後，故障的案例相對而言還是太少，是以故障率的不确定性還非常高。

當一個組織或公司向把離線資料架構轉換為線上流式處理模式，他們應該已經積累的不少曆史資料。這隻是關于故障率的先驗信念，可以作為故障率的先驗分布注入。

伽馬分布可以友善的和泊松分布

共轭(conjugate)

，是以我們可以使用伽馬分布去代表先驗機率，并且在觀測到一些資料之後，可以很簡答的計算得到後驗伽馬分布（不需要使用很複雜的數學計算）。因為伽馬分布的 ɑ 和 β 兩個參數可以分别被解釋為故障數目和運作天數，是以這些值可以通過曆史資料計算出來，并通過Spark Structured Streaming作業持續的計算和更新。請看示例：

priors_df = spark.read.\
  option("header", True).\
  option("inferSchema", True).\
  csv(data_path + "/historical/").\
  select("model", "failure").\
  groupBy("model").\
  agg({"*": "count", "failure": "avg"}).\
  withColumn("alpha", col("count(1)") * col("avg(failure)")).\
  withColumn("beta", col("count(1)")).\
  select("model", "alpha", "beta").\
  cache()

agg_stream = spark.readStream.\
  schema(raw_df.schema).\
  option("header", True).\
  option("maxFilesPerTrigger", 1).\
  csv(data_path + "/current/").\
  select("model", "failure").\
  groupBy("model").\
  agg({"*": "count", "failure": "avg"}).\
  filter("avg(failure) > 0").\
  orderBy(desc("avg(failure)")).\
  limit(6)

max_avg_failure = priors_df.select(avg(col("alpha") /
col("beta")).alias("failure_rate")).agg(max("failure_rate")).collect()[0][0]

def to_gamma_pdf(count, avgf, prior_alpha, prior_beta):
  # Note: priors are added below
  k = prior_alpha + count * avgf
  theta = 1.0 / (prior_beta + count)
  rates = np.geomspace(max_avg_failure / 1000, max_avg_failure, 100)
  pdfs = stats.gamma.pdf(rates, a=k, scale=theta)
  return list(zip(rates.tolist(), pdfs.tolist()))

to_gamma_pdf_udf = udf(to_gamma_pdf, to_gamma_schema)

gamma_stream = agg_stream.\
  join(priors_df, "model", how='left').\
  fillna(0, subset=["alpha", "beta"]).\
  withColumn("gamma", to_gamma_pdf_udf("count(1)", "avg(failure)", "alpha", "beta")).\
  select("model", explode("gamma").alias("rates_pdfs")).\
  select("model", col("rates_pdfs.rates").alias("rate"), col("rates_pdfs.pdfs").alias("pdf"))           

通過曆史資料計算産生了各個硬碟型号先驗分布的ɑ和β參數，然後可以很友善的和流式資料進行Join（這裡的流式資料是我們模拟得到的）

開始運作之後，程式生成了一個不斷更新的圖表來展示目前最高故障率的硬碟型号的故障率分布情況。它可以很容易地生成報警，或者根據機率分布的其他功能生成其他警報，并且可以使用幾乎相同的代碼生成簡單的離線分析程式。

結語

當你在拟合一個分布曲線時，實際上發生的事情比你想象的更多。你可能已經在不知不覺中使用了“機率派”或“貝葉斯派”的推理方法。貝葉斯派觀點有一些優點，因為在尋找答案過程中考慮了先驗知識，并産生的是答案的一個機率分布而不隻是一個估計值。通過Spark的Structured Streaming技術以及對先驗和後驗分布的了解，可以将貝葉斯方法論友善的應用到現實世界的實踐中。