贝叶斯方法应用：检测时间序列拐点

随着时间推移，制造设备比如贴片机的位置由于各种原因会产生小的偏差。这些偏差可能是阶跃，也有可能是渐变的形式。由于偏差值很小，产线的自动光学检测设备并不会报警；然而小的偏差如果不经处理，经过一定时间累积会产生较大偏差，影响产品质量。

为了能够提前发现设备偏差并在产生质量问题以前及时调整，产线工程师需要有工具能够持续监测设备各个部件的性能，在出现平均值或标准偏差的拐点时提醒工程师及时做设备调整。

这是一类时间序列检测拐点的问题，可以通过贝叶斯方法求解任意时间点出现拐点的概率。

贝叶斯方法基础

摘抄自刘未鹏：平凡而又神奇的贝叶斯方法

贝叶斯公式是怎么来的？

使用 wikipedia 上的一个例子：

一所学校里面有 60% 的男生，40% 的女生。男生总是穿长裤，女生则一半穿长裤一半穿裙子。有了这些信息之后我们可以容易地计算“随机选取一个学生，他（她）穿长裤的概率和穿裙子的概率是多大”，这个就是前面说的“正向概率”的计算。然而，假设你走在校园中，迎面走来一个穿长裤的学生（很不幸的是你高度近似，你只看得见他（她）穿的是否长裤，而无法确定他（她）的性别），你能够推断出他（她）是男生的概率是多大吗？

一些认知科学的研究表明（《决策与判断》以及《Rationality for Mortals》第12章：小孩也可以解决贝叶斯问题），我们对形式化的贝叶斯问题不擅长，但对于以频率形式呈现的等价问题却很擅长。在这里，我们不妨把问题重新叙述成：你在校园里面随机游走，遇到了 N 个穿长裤的人（仍然假设你无法直接观察到他们的性别），问这 N 个人里面有多少个女生多少个男生。

假设学校里面人的总数是 U 个。60% 的男生都穿长裤，于是我们得到了 U P(Boy) P(Pants|Boy) 个穿长裤的（男生）（其中 P(Boy) 是男生的概率 = 60%，这里可以简单的理解为男生的比例；P(Pants|Boy) 是条件概率，即在 Boy 这个条件下穿长裤的概率是多大，这里是 100% ，因为所有男生都穿长裤）。40% 的女生里面又有一半（50%）是穿长裤的，于是我们又得到了 U P(Girl) P(Pants|Girl) 个穿长裤的（女生）。加起来一共是 U P(Boy) P(Pants|Boy) + U P(Girl) P(Pants|Girl) 个穿长裤的，其中有 U P(Girl) P(Pants|Girl) 个女生。两者一比就是你要求的答案。

下面我们把这个答案形式化一下：我们要求的是 P(Girl|Pants) （穿长裤的人里面有多少女生），我们计算的结果是 U P(Girl) P(Pants|Girl) / [U P(Boy) P(Pants|Boy) + U P(Girl) P(Pants|Girl)] 。容易发现这里校园内人的总数是无关的，可以消去。于是得到

P(Girl|Pants) = P(Girl) P(Pants|Girl) / [P(Boy) P(Pants|Boy) + P(Girl) * P(Pants|Girl)]

注意，如果把上式收缩起来，分母其实就是 P(Pants) ，分子其实就是 P(Pants, Girl) 。而这个比例很自然地就读作：在穿长裤的人（ P(Pants) ）里面有多少（穿长裤）的女孩（ P(Pants, Girl) ）。

上式中的 Pants 和 Boy/Girl 可以指代一切东西，所以其一般形式就是：

$$P(B|A)=P(A|B)P(B)/[P(A|B)P(B)+P(A|B)\ast P(B)]$$

收缩起来就是 $P(B|A)=P(AB)/P(A)$

其实这个就等于 $P(B|A)\ast P(A)=P(AB)$

产线设备性能拐点案例：数学建模

以贴片机为例：根据业务场景，确定了三类需要检测的性能拐点：

• 贴片位置偏差平均值阶跃变化

• 贴片位置偏差平均值渐变

• 贴片位置浮动大小（标准偏差）突变

根据这三类性能拐点，将问题描述为包含拐点的混合线性模型：

$$y(t)={\beta }_0{\varphi }_{-}^{\theta }+{\beta }_1{\zeta }_{-}^{\theta }+{\beta }_2{\zeta }_{+}^{\theta }+{\beta }_3{\varphi }_{+}^{\theta }+\xi (t)$$

其中 $\varphi$ 和 $\zeta$都是阶跃函数，$\varphi$为常数阶跃，$\zeta$为线性阶跃：

$${\varphi }_{-}^{\theta }=\{\begin{array}{ll}1& \text{if }t\text{ <= }\theta \\ 0& \text{else}\end{array}$$

$${\varphi }_{+}^{\theta }=\{\begin{array}{ll}1& \text{if }t\text{ >= }\theta \\ 0& \text{else}\end{array}$$

$${\zeta }_{-}^{\theta }=\{\begin{array}{ll}\theta -t& \text{if }t\text{ <= }\theta \\ 0& \text{else}\end{array}$$

$${\zeta }_{+}^{\theta }=\{\begin{array}{ll}\theta -t& \text{if }t\text{ >= }\theta \\ 0& \text{else}\end{array}$$

计算${\varphi }_{-}^{\theta }$，${\varphi }_{+}^{\theta }$，${\zeta }_{-}^{\theta }$，${\zeta }_{+}^{\theta }$的Python代码示例如下：

import numpy as np
def xiMinus(theta, t, mode = 'constant'):
    scale = 1
    if mode == 'constant':
        if t <= theta:
            return -1.0 / scale
        else:
            return 0.0
    if mode == 'linear':
        if t <= theta:
            return 1.0 / scale * (theta - t)
        else:
            return 0.0

def xiPlus(theta, t, mode = 'constant'):
    scale = 1
    if mode == 'constant':
        if t <= theta:
            return 0.0
        else:
            return 1.0 / scale
    if mode == 'linear':
        if t <= theta:
            return 0.0
        else:
            return 1.0 / scale * (t - theta)

$\xi (t)$ 为在平均值左右正态分布的随机浮动，分布的标准偏差也在拐点测试的范围内，i.e.假设数据浮动的幅度也有可能出现阶跃。所以定义随机浮动的标准偏差：

$$STD(\xi (t))=\sigma (1+$$