此页面由 Cloud Translation API 翻译。
Switch to English

多变化点检测和贝叶斯模型选择

贝氏模型选择

在TensorFlow.org上查看 在Google Colab中运行 在GitHub上查看源代码 下载笔记本

进口货

import numpy as np
import tensorflow.compat.v2 as tf
tf.enable_v2_behavior()
import tensorflow_probability as tfp
from tensorflow_probability import distributions as tfd

from matplotlib import pylab as plt
%matplotlib inline
import scipy.stats

任务:具有多个变更点的变更点检测

考虑一个变更点检测任务:事件发生的速率随着时间的变化而变化,这是由某些生成数据的系统或进程的(未观察到的)状态的突然变化驱动的。

例如,我们可能会观察到一系列如下计数:

true_rates = [40, 3, 20, 50]
true_durations = [10, 20, 5, 35]

observed_counts = np.concatenate([
  scipy.stats.poisson(rate).rvs(num_steps)
    for (rate, num_steps) in zip(true_rates, true_durations)
]).astype(np.float32)

plt.plot(observed_counts)
[<matplotlib.lines.Line2D at 0x7fe01d3c2390>]

png

这些可能代表数据中心中的故障数量,网页的访问者数量,网络链接上的数据包数量等。

请注意,仅查看数据并不清楚有多少不同的系统状态。您能否确定三个切换点中的每一个出现在哪里?

已知状态数

我们将首先考虑(可能不现实)先验未知状态数的情况。在这里,我们假设我们知道有四个潜在状态。

我们将此问题建模为一个切换(非均匀的)泊松过程:在每个时间点,发生的事件数都是泊松分布的,并且事件的发生率由未观察到的系统状态$ z_t $决定:

$$x_t \sim \text{Poisson}(\lambda_{z_t})$$

潜在状态是离散的:$ z_t \ in {1,2,3,4} $,因此$ \ lambda = [\ lambda_1,\ lambda_2,\ lambda_3,\ lambda_4] $是一个简单的向量,每个向量都包含泊松率州。为了模拟状态随时间的演变,我们将定义一个简单的过渡模型$ p(z_t | z_ {t-1})$:假设在每一步中,我们都以一定的概率$ p $停留在先前状态,并以概率$ 1-p $随机地均匀地过渡到另一个状态。初始状态也是随机选择的,因此我们有:

$$ \begin{align*} z_1 &\sim \text{Categorical}\left(\left\{\frac{1}{4}, \frac{1}{4}, \frac{1}{4}, \frac{1}{4}\right\}\right)\\ z_t | z_{t-1} &\sim \text{Categorical}\left(\left\{\begin{array}{cc}p & \text{if } z_t = z_{t-1} \\ \frac{1-p}{4-1} & \text{otherwise}\end{array}\right\}\right) \end{align*}$$

这些假设对应于具有泊松排放的隐马尔可夫模型 。我们可以使用tfd.HiddenMarkovModel在TFP中对它们进行编码。首先,我们定义转换矩阵和初始状态的统一先验:

num_states = 4

initial_state_logits = np.zeros([num_states], dtype=np.float32) # uniform distribution

daily_change_prob = 0.05
transition_probs = daily_change_prob / (num_states-1) * np.ones(
    [num_states, num_states], dtype=np.float32)
np.fill_diagonal(transition_probs,
                 1-daily_change_prob)

print("Initial state logits:\n{}".format(initial_state_logits))
print("Transition matrix:\n{}".format(transition_probs))
Initial state logits:
[0. 0. 0. 0.]
Transition matrix:
[[0.95       0.01666667 0.01666667 0.01666667]
 [0.01666667 0.95       0.01666667 0.01666667]
 [0.01666667 0.01666667 0.95       0.01666667]
 [0.01666667 0.01666667 0.01666667 0.95      ]]

接下来,我们使用可训练变量表示与每个系统状态关联的速率,构建tfd.HiddenMarkovModel分布。我们在对数空间中对比率进行参数化,以确保它们为正值。

# Define variable to represent the unknown log rates.
trainable_log_rates = tf.Variable(
  np.log(np.mean(observed_counts)) + tf.random.normal([num_states]),
  name='log_rates')

hmm = tfd.HiddenMarkovModel(
  initial_distribution=tfd.Categorical(
      logits=initial_state_logits),
  transition_distribution=tfd.Categorical(probs=transition_probs),
  observation_distribution=tfd.Poisson(log_rate=trainable_log_rates),
  num_steps=len(observed_counts))

最后,我们定义模型的总对数密度,包括在速率之前先提供弱信息的LogNormal,然后运行优化程序以计算与观察到的计数数据拟合的最大后验 (MAP)。

rate_prior = tfd.LogNormal(5, 5)

def log_prob():
 return (tf.reduce_sum(rate_prior.log_prob(tf.math.exp(trainable_log_rates))) +
         hmm.log_prob(observed_counts))

optimizer = tf.keras.optimizers.Adam(learning_rate=0.1)

@tf.function(autograph=False)
def train_op():
  with tf.GradientTape() as tape:
    neg_log_prob = -log_prob()
  grads = tape.gradient(neg_log_prob, [trainable_log_rates])[0]
  optimizer.apply_gradients([(grads, trainable_log_rates)])
  return neg_log_prob, tf.math.exp(trainable_log_rates)
for step in range(201):
  loss, rates = [t.numpy() for t in train_op()]
  if step % 20 == 0:
    print("step {}: log prob {} rates {}".format(step, -loss, rates))

print("Inferred rates: {}".format(rates))
print("True rates: {}".format(true_rates))
step 0: log prob -539.8363037109375 rates [ 60.120716  46.096066 115.06791   21.380898]
step 20: log prob -251.49668884277344 rates [49.44194  22.835463 42.680107  4.509814]
step 40: log prob -244.906982421875 rates [49.37409   22.341614  39.15951    2.7087321]
step 60: log prob -244.69183349609375 rates [49.476925  21.316504  40.08259    2.7874506]
step 80: log prob -244.27362060546875 rates [49.519707  22.292286  40.211983   3.0713978]
step 100: log prob -244.25796508789062 rates [49.507477  22.055244  40.16306    3.1456223]
step 120: log prob -244.2534942626953 rates [49.510178  21.981022  40.136353   3.1146054]
step 140: log prob -244.2531280517578 rates [49.514404  22.00701   40.136944   3.1029546]
step 160: log prob -244.25303649902344 rates [49.513     22.018309  40.144352   3.1072166]
step 180: log prob -244.25303649902344 rates [49.51299   22.01862   40.142887   3.1083844]
step 200: log prob -244.2530517578125 rates [49.51313   22.017632  40.142963   3.1076818]
Inferred rates: [49.51313   22.017632  40.142963   3.1076818]
True rates: [40, 3, 20, 50]

有效!请注意,此模型中的潜在状​​态仅在排列之前是可识别的,因此我们恢复的速率处于不同的顺序,并且有一些噪音,但是通常它们匹配得很好。

恢复状态轨迹

现在,我们已经适应的模式,我们可能要重新构建状态模型认为,该系统是在每个时间步长。

这是一个后验任务:给定观察到的计数$ x_ {1:T} $和模型参数(比率)$ \ lambda $,我们要根据后验分布$ p(z_ { 1:T} | x_ {1:T},\ lambda)$。在隐马尔可夫模型中,我们可以使用标准消息传递算法有效地计算此分布的边际和其他属性。特别地, posterior_marginals方法将在每个时间步长$ t $上有效计算(使用前向后退算法 )离散离散状态$ Z_t $上的边际概率分布$ p(Z_t = z_t | x_ {1:T})$ 。

# Runs forward-backward algorithm to compute marginal posteriors.
posterior_dists = hmm.posterior_marginals(observed_counts)
posterior_probs = posterior_dists.probs_parameter().numpy()

通过绘制后验概率,我们可以恢复模型对数据的“解释”:每个状态在哪个时间点处于活动状态?

03

png

在这种(简单)情况下,我们看到该模型通常是非常有信心的:在大多数时间步中,它基本上将所有概率质量分配给四个状态中的一个。幸运的是,这些解释看起来很合理!

我们还可以根据每个时间步上与最可能的潜在状态相关的速率来可视化此后验,将概率后验简化为一个单一的解释:

most_probable_states = np.argmax(posterior_probs, axis=1)
most_probable_rates = rates[most_probable_states]
fig = plt.figure(figsize=(10, 4))
ax = fig.add_subplot(1, 1, 1)
ax.plot(most_probable_rates, c='green', lw=3, label='inferred rate')
ax.plot(observed_counts, c='black', alpha=0.3, label='observed counts')
ax.set_ylabel("latent rate")
ax.set_xlabel("time")
ax.set_title("Inferred latent rate over time")
ax.legend(loc=4)
<matplotlib.legend.Legend at 0x7fe01d0afda0>

png

技术说明:不是每个单独的时间步长最可能的状态,而是$ z ^ * _ t = \ text {argmax} _ {z_t} p(z_t | x_ {1:T})$,我们可能会要求最可能的状态潜在轨迹 ,$ z ^ * = \ text {argmax} _z p(z | x_ {1:T})$(甚至包括来自后轨迹的样本!),同时考虑了时间步长之间的相关性。为了说明这种差异,假设剪刀石头布的玩家在40%的时间里都玩石头,但从来没有连续两次玩过石头:在每个时间点上,石头都是最可能的边缘状态,但是“石头,石头,石头”。绝对不是最可能的轨迹-实际上,它的可能性为零!

TODO(davmre):一旦tfp.HiddenMarkovModel实现了Viterbi算法以找到最高概率的轨迹,请更新此部分以使用它。

状态数未知

在实际问题中,我们可能不知道正在建模的系统中状态的“真实”数目。这可能并不总是一个问题:如果您不特别关心未知状态的标识,则可以运行一个模型,该模型具有比您知道模型所需的状态更多的状态,并学习(类似)一堆重复项实际状态的副本。但让我们假设您确实关心推断潜在状态的“真实”数目。

我们可以将其视为贝叶斯模型选择的一种情况:我们有一组候选模型,每个候选模型具有不同数量的潜在状态,并且我们想选择最有可能生成观测数据的模型。为此,我们计算了每个模型下数据的边际可能性(我们也可以在模型本身上添加先验,但是在此分析中这不是必需的; 贝叶斯Occam的剃刀证明足以编码一个偏爱简单的模型)。

不幸的是,真实的边际可能性在离散状态$ z_ {1:T} $和速率参数(\ vector)的两个向量上积分,

$$p(x_{1:T}) = \int p(x_{1:T}, z_{1:T}, \lambda) dz d\lambda,$$

对于此模型而言是不易处理的。为了方便起见,我们将使用所谓的“ 经验贝叶斯 ”或“ II型最大似然”估计值对其进行近似:我们将完全整合与每个系统状态相关的(未知)速率参数$ \ lambda $,而不是优化其价值:

$$\tilde{p}(x_{1:T}) = \max_\lambda \int p(x_{1:T}, z_{1:T}, \lambda) dz$$

这种近似可能会过拟合,即,与真正的边际可能性相比,它更喜欢更复杂的模型。我们可以考虑更忠实的近似值,例如,优化变分下界,或使用蒙特卡洛估计器(例如, 退火的重要性采样) ;这些(严重)超出了本笔记本的范围。 (有关贝叶斯模型选择和近似的更多信息,优秀的机器学习:概率观点的第7章是很好的参考。)

原则上,我们可以简单地通过多次使用不同num_states值重新运行优化来进行模型比较,但这将需要大量工作。在这里,我们将展示如何使用TFP的batch_shape机制进行矢量化并行考虑多个模型。

过渡矩阵和初始状态优先级:现在,我们将构建一批过渡矩阵和优先级逻辑,而不是构建单个模型描述,对于每个候选模型(最大max_num_states 。为了便于批处理,我们需要确保所有计算都具有相同的“形状”:这必须与我们适合的最大模型的尺寸相对应。为了处理较小的模型,我们可以将它们的描述“嵌入”状态空间的最高维度,从而将其余维度有效地视为从未使用过的伪状态。

max_num_states = 10

def build_latent_state(num_states, max_num_states, daily_change_prob=0.05):

  # Give probability exp(-100) ~= 0 to states outside of the current model.
  initial_state_logits = -100. * np.ones([max_num_states], dtype=np.float32)
  initial_state_logits[:num_states] = 0.

  # Build a transition matrix that transitions only within the current
  # `num_states` states.
  transition_probs = np.eye(max_num_states, dtype=np.float32)
  if num_states > 1:
    transition_probs[:num_states, :num_states] = (
        daily_change_prob / (num_states-1))
    np.fill_diagonal(transition_probs[:num_states, :num_states],
                     1-daily_change_prob)
  return initial_state_logits, transition_probs

# For each candidate model, build the initial state prior and transition matrix.
batch_initial_state_logits = []
batch_transition_probs = []
for num_states in range(1, max_num_states+1):
  initial_state_logits, transition_probs = build_latent_state(
      num_states=num_states,
      max_num_states=max_num_states)
  batch_initial_state_logits.append(initial_state_logits)
  batch_transition_probs.append(transition_probs)

batch_initial_state_logits = np.array(batch_initial_state_logits)
batch_transition_probs = np.array(batch_transition_probs)
print("Shape of initial_state_logits: {}".format(batch_initial_state_logits.shape))
print("Shape of transition probs: {}".format(batch_transition_probs.shape))
print("Example initial state logits for num_states==3:\n{}".format(batch_initial_state_logits[2, :]))
print("Example transition_probs for num_states==3:\n{}".format(batch_transition_probs[2, :, :]))
Shape of initial_state_logits: (10, 10)
Shape of transition probs: (10, 10, 10)
Example initial state logits for num_states==3:
[   0.    0.    0. -100. -100. -100. -100. -100. -100. -100.]
Example transition_probs for num_states==3:
[[0.95  0.025 0.025 0.    0.    0.    0.    0.    0.    0.   ]
 [0.025 0.95  0.025 0.    0.    0.    0.    0.    0.    0.   ]
 [0.025 0.025 0.95  0.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    1.    0.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    1.    0.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    1.    0.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    1.    0.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    1.    0.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    1.    0.   ]
 [0.    0.    0.    0.    0.    0.    0.    0.    0.    1.   ]]

现在,我们与上述类似地进行。这次,我们将在trainable_rates使用一个额外的批处理维度,以分别适合所考虑的每个模型的费率。

trainable_log_rates = tf.Variable(
    (np.log(np.mean(observed_counts)) *
     np.ones([batch_initial_state_logits.shape[0], max_num_states]) +
     tf.random.normal([1, max_num_states])),
     name='log_rates')
    
hmm = tfd.HiddenMarkovModel(
  initial_distribution=tfd.Categorical(
      logits=batch_initial_state_logits),
  transition_distribution=tfd.Categorical(probs=batch_transition_probs),
  observation_distribution=tfd.Poisson(log_rate=trainable_log_rates),
  num_steps=len(observed_counts))

在计算总对数概率时,我们小心地仅总结每个模型组件实际使用的费率的先验值:

rate_prior = tfd.LogNormal(5, 5)

def log_prob():
  prior_lps = rate_prior.log_prob(tf.math.exp(trainable_log_rates))
  prior_lp = tf.stack(
      [tf.reduce_sum(prior_lps[i, :i+1]) for i in range(max_num_states)])
  return prior_lp + hmm.log_prob(observed_counts)
@tf.function(autograph=False)
def train_op():
  with tf.GradientTape() as tape:
    neg_log_prob = -log_prob()
  grads = tape.gradient(neg_log_prob, [trainable_log_rates])[0]
  optimizer.apply_gradients([(grads, trainable_log_rates)])
  return neg_log_prob, tf.math.exp(trainable_log_rates)

现在,我们优化已构建的批次目标,同时拟合所有候选模型:

for step in range(201):
  loss, rates =  [t.numpy() for t in train_op()]
  if step % 20 == 0:
    print("step {}: loss {}".format(step, loss))
step 0: loss [932.1646  832.1462  828.4412  837.84705 846.1599  417.9551  419.0383
 417.9084  425.87964 366.82443]
step 20: loss [807.2557  286.25125 257.90186 256.22498 262.71576 270.11447 279.765
 267.6623  276.87903 284.21024]
step 40: loss [803.20874 279.9999  254.07617 254.04176 262.7116  257.1419  265.64944
 263.69067 270.1178  278.18292]
step 60: loss [803.0721  271.3898  244.30405 244.61728 252.61594 256.8096  263.2815
 261.81732 268.69806 277.49   ]
step 80: loss [802.9     271.44714 244.07954 244.36232 252.20078 256.4545  261.7729
 260.18414 267.57748 276.11465]
step 100: loss [802.9022  271.3092  243.98022 244.25925 251.41597 256.0926  260.91467
 257.74948 265.95273 275.0009 ]
step 120: loss [802.8999  271.30957 243.9751  244.25511 250.79546 255.72038 260.2765
 256.68982 264.83832 274.82697]
step 140: loss [802.89966 271.3077  243.97357 244.25308 250.65877 255.32945 259.61533
 255.69922 264.014   274.46405]
step 160: loss [802.8996  271.3077  243.97354 244.25311 250.64474 254.89218 258.482
 254.8172  263.24512 273.9937 ]
step 180: loss [802.8997  271.3077  243.97356 244.25302 250.64305 253.76013 257.86176
 254.0182  262.5052  273.49908]
step 200: loss [802.89954 271.3077  243.9734  244.2529  250.64278 253.11316 257.45514
 253.2917  261.7873  273.0126 ]

num_states = np.arange(1, max_num_states+1)
plt.plot(num_states, -loss)
plt.ylim([-400, -200])
plt.ylabel("marginal likelihood $\\tilde{p}(x)$")
plt.xlabel("number of latent states")
plt.title("Model selection on latent states")
Text(0.5, 1.0, 'Model selection on latent states')

png

检查可能性,我们发现(近似)边际可能性更喜欢三态或四态模型(具体顺序可能在此笔记本的运行之间有所不同)。这似乎很合理-“真实”模型具有四个状态,但是仅查看数据就很难排除三状态的解释。

我们还可以提取适合每个候选模型的费率:

for i, learned_model_rates in enumerate(rates):
  print("rates for {}-state model: {}".format(i+1, learned_model_rates[:i+1]))
rates for 1-state model: [32.98682]
rates for 2-state model: [44.949432   3.1397576]
rates for 3-state model: [22.019554   3.1073658 47.4462   ]
rates for 4-state model: [22.01749    3.1073527 49.513165  40.14324  ]
rates for 5-state model: [22.016987   3.1073353 49.470993  40.127758  49.488205 ]
rates for 6-state model: [22.016418    0.12033073 50.25381    40.12794    48.88818     3.1076846 ]
rates for 7-state model: [4.9506187e+01 2.2016148e+01 4.9468941e+01 4.7518797e-03 4.9311584e+01
 3.1077573e+00 4.0113823e+01]
rates for 8-state model: [4.0115150e+01 4.3629836e-02 4.9482445e+01 4.5004072e-05 3.1080871e+00
 5.0322604e-01 4.9483521e+01 2.2015779e+01]
rates for 9-state model: [4.0034302e+01 7.8987077e-02 4.9487354e+01 3.3131179e-03 4.0034004e+01
 4.7514364e-01 4.9488628e+01 2.2016052e+01 3.1080682e+00]
rates for 10-state model: [39.950623    3.0235524  39.950375    0.16000797 39.950424   21.830935
 21.831202    3.0232046  49.51654     3.0235553 ]

并绘制每个模型为数据提供的解释:

posterior_probs = hmm.posterior_marginals(
    observed_counts).probs_parameter().numpy()
most_probable_states = np.argmax(posterior_probs, axis=-1)
fig = plt.figure(figsize=(14, 12))
for i, learned_model_rates in enumerate(rates):
  ax = fig.add_subplot(4, 3, i+1)
  ax.plot(learned_model_rates[most_probable_states[i]], c='green', lw=3, label='inferred rate')
  ax.plot(observed_counts, c='black', alpha=0.3, label='observed counts')
  ax.set_ylabel("latent rate")
  ax.set_xlabel("time")
  ax.set_title("{}-state model".format(i+1))
  ax.legend(loc=4)
plt.tight_layout()

png

很容易看出一态,二态和(更巧妙)的三态模型提供的解释不足。有趣的是,以上四个状态的所有模型都提供了基本相同的解释!这可能是因为我们的“数据”相对干净,几乎没有余地可供选择。在更现实的数据上,我们希望容量更高的模型能够逐步更好地拟合数据,并在某些折衷点上,模型的复杂性会超出改进的拟合度。

扩展名

该笔记本中的模型可以通过多种方式直接扩展。例如:

  • 允许潜在状态具有不同的概率(某些状态可能很常见,也很罕见)
  • 允许在潜在状态之间进行非均匀转换(例如,了解到机器崩溃通常伴随着系统重新引导,通常伴随着一段良好的性能,等等)
  • 其他排放模型,例如NegativeBinomial用于对计数数据中的变化离散进行建模,或者对实值数据进行连续分布,例如Normal