学习笔记 | 贝叶斯块算法

施工中

获取探测器的计数随时间或其他物理量的变化可以帮助了解发生的物理过程的各种性质。例如在 X 射线和 γ 射线天文学中十分关注光子计数随时间的变化，这一类变化被称为光变曲线，可以借此计算相应物理过程的空间尺度等性质。

基础知识¶

核辐射测量的统计性质¶

探测 X/γ 射线需要使用核辐射探测器，每探测到一个事例，探测器将会记录相应事例的物理信息，如时间、能量、径迹等。核辐射的测量充满了随机性，但是可以用统计分布来描述其中的随机性。单个放射性粒子的衰变过程是一个伯努利事件，其衰变常数为 $\lambda$，则在时间 $0\sim t$ 内发生衰变的可能性为：

\begin{align} p = 1 - e^{-\lambda t} \end{align}

对于由 $n_0$ 个放射性粒子组成的体系，彼此之间发生衰变是独立的，则体系的衰变过程是一个 $n_0$ 重伯努利过程，即 $t$ 时刻发生了衰变的粒子数目 $n$ 满足二项分布：

\begin{equation}\begin{split} P(n|n_0) &= C_{n_0}^np^n(1-p)^{n_0-n} \\ \mu &= E(n) = n_0p \\ \sigma^2 &= D(n) = n_0p(1-p) \end{split}\end{equation}

当 $n_0$ 足够大，而 $p$ 足够小时，二项分布可以近似为泊松分布：

\begin{equation}\begin{split} P(n) &= \frac{\mu^n}{n!}e^{-\mu} \\ \mu &= n_0p \\ \sigma &= \sqrt{\mu} \end{split}\end{equation}

对于一般的物理过程而言，都满足该条件。当均值 $\mu$ 足够大时，泊松分布可以进一步近似为高斯分布。在辐射测量中，我们获取到的是探测时间 $\mathrm{d}t$ 内的总计数，假设 $\mathrm{d}t$ 时间内源的强度 $\lambda$ （单位为 $\mathrm{s}^{-1}$）不变，则有：

$$\mu = \lambda\mathrm{d}t$$

两个相邻事例的时间间隔为 $T$ 表示在第一个事例后 $0\sim T$ 时间内没有事例，而在 $T\sim T+\mathrm{d}t$ 时间内有一个事例，其概率可以表示为：

\begin{equation}\begin{split} P(t\le T \le t+\mathrm{d}t) &= P_t(0)\times P_{\mathrm{d}t}(1)\\ P_t(0) = \frac{(\lambda t)^0}{0!}e^{-\lambda t} &= e^{-\lambda t} \\ P_{\mathrm{d}t}(1) = \frac{(\lambda\mathrm{d}t)^1}{1!}e^{-\lambda\mathrm{d}t} &\approx \lambda\mathrm{d}t\quad(e^{-\lambda\mathrm{d}t} \to 1) \end{split}\end{equation}

设 $f(t)$ 为概率密度函数，则有：

\begin{align} f(t)\mathrm{d}t = P(t\le T \le t+\mathrm{d}t) = \lambda e^{-\lambda t}\mathrm{d}t \end{align}

可见时间间隔 $T$ 满足指数分布，可以求得其期望值 $\overline{t}=\frac{1}{\lambda}$。

贝叶斯定理¶

对于一个模型 $\mathscr{M}$，假设其所有参数构成向量 $\theta$。则根据已知数据 $D$，贝叶斯定理可以表示为：

\begin{align} P(\theta|D,\mathscr{M}) = \frac{P(D|\theta,\mathscr{M})P(\theta|\mathscr{M})}{P(D|\mathscr{M})} \end{align}

eq-bayesian

其中：

• $P(\theta|D,\mathscr{M})$ 为 $\theta$ 的后验概率（已知 $D$ 情况下 $\theta$ 的条件概率）
• $P(\theta|\mathscr{M})$ 为 $\theta$ 的先验概率
• $P(D|\theta,\mathscr{M})$ 为 $D$ 的后验概率或者在特定 $D$ 下，$\theta$ 的似然性，即 $L(\theta|D,\mathscr{M})$
• $P(D|\mathscr{M})$ 为 $D$ 的先验概率

对于一组模型 $\mathscr{M}_1,\mathscr{M}_2,\dots,\mathscr{M}_n$，所有的约束、相关背景知识或者假设被称为信息背景，记作 $I$。则在给定 $D$ 和 $I$ 的情况下，不同模型的后验概率可以表示为：

\begin{align} P(\mathscr{M}_k|D,I) = \frac{P(D|\mathscr{M}_k,I)P(\mathscr{M}_k|I)}{P(D|I)} \end{align}

信息背景应当是不变的，所以我们可以省略。因此两种模型对数据的符合程度可以使用其后验概率的比值（即 odds ratio）表示：

\begin{align} \mathscr{O}_{ij} = \frac{P(\mathscr{M}_i|D)}{P(\mathscr{M}_j|D)} = \frac{P(D|\mathscr{M}_i)P(\mathscr{M}_i)}{P(D|\mathscr{M}_j)P(\mathscr{M}_j)} \end{align}

根据贝叶斯公式，由于 $P(\theta_k|D,\mathscr{M}_k)$ 归一化，由式 (eq-bayesian) 可以得到：

\begin{align} \mathscr{L}(\mathscr{M}_k|D) \equiv P(D|\mathscr{M}_k) = \int P(D|\theta_k,\mathscr{M}_k)P(\theta_k|\mathscr{M}_k)\mathrm{d}\theta_k \\ J(\mathscr{M}_k|D) \equiv P(\mathscr{M}_k)\int P(D|\theta_k,\mathscr{M}_k)P(\theta_k|\mathscr{M}_k)\mathrm{d}\theta_k \end{align}

eq-global-likelihood

其中 $\mathscr{L}(\mathscr{M}_k|D)$ 被称为全局似然函数（global likelihood，叫边缘概率更多，marginal probability），$J(\mathscr{M}_k|D)$ 是数据和模型的联合概率（joint probability），包含了先验信息而与模型参数无关，则有：

\begin{align} \mathscr{O}_{ij} = \frac{\mathscr{L}(\mathscr{M}_i|D)}{\mathscr{L}(\mathscr{M}_j|D)}\rho \end{align}

eq-odds-ratio

其中：

\begin{align} \rho = \frac{P(\mathscr{M}_i)}{P(\mathscr{M}_j)} \end{align}

光变曲线¶

一定时间范围内的计数随时间的变化曲线就是我们所需的光变曲线，因此光变曲线可以理解为探测事例根据时间绘制的频数分布直方图。直方图的每一组在这里被称为分箱（bin），组距被称为 bin 宽，各组的界限值被称为 bin edge，直方图一个很重要的特性的分箱之间是等距的。

对事例进行分箱的会丢弃大量信息，而且分箱的大小和位置对最终分析结果的影响很大。bin 宽越小，能够展现更精细的结构，但是单一分箱内的计数过少，统计性不足；bin 宽越大，统计性越强，但是丢失的物理信息更多。所以我们需要一种方法，既能够保证统计性，同时还能尽可能地展现精细的结构。这种方法应当允许分箱拥有不同的 bin 宽，而且 bin edge 能够根据数据自适应确定。

一个简单的想法是遍历每一种可能的分箱方案，而可能方案的总数为 $2^N$ ，其中 $N$ 为数据点的数目。显而易见随着数据量的增大，这种方法很快变得不可用。Scargle 在 1998 年提出了基于似然函数的迭代算法⁽¹⁾，能够在可接受的时间范围内得出近似解；而在 2013 年，Scargle 基于动态规划（Dynamic Programming）和块评价函数提出了一种新算法⁽²⁾，能够在 $O(N^2)$ 的时间内求出最优解。在本文中，前者称为近似算法，后者为DP算法，两者统称为贝叶斯块算法。

贝叶斯块算法¶

贝叶斯块算法将光变曲线表示为亮度随时间的分段常数，可以用于估计背景水平以及物理事件的时间尺度、空间尺度与强度等信息。

在该算法中，单个的数据点被称为 Cell，而对数据进行分割的点被称为分割点 Change Point（简写为 cp），而相邻两个分割点之间（包括第一个分割点之前与最后一个分割点之后的）的所有数据点的集合被称为块 Block，某种特定的分割方案（或分割点的组合）被称为组合 Partition。

算法的关键思想是不同块之间的评价是独立的，仅取决于该块的性质（例如块的长度、数据的情况等）。例如简单的分段常数模型具有这些参数：

• $N_\mathrm{cp}$：分割点的数目
• $t^\mathrm{cp}_k$：第 $k$ 个块的起始点（即第 $k-1$ 个和第 $k$ 个块的分割点）
• $X_k$：第 $k$ 个块的信号强度

其中 $k\in\{1,2,\dots,N_\mathrm{cp}+1\}$，块最重要的两个参数是强度（即 $X_k$）与长度（起始点与下一个分割点之间的间隔）。模型关键的问题在于要分多少个块，近似算法和DP算法对此的处理方式有所区别，具体参见后文。

数据模式¶

Time-Tagged Event Data¶

Time-Tagged Event（时间标记事件，TTE）是指记录事例入射时间的模式，表示为：

\begin{align} D_\mathrm{TTE}: \{t_n,n=1\dots N\} \end{align}

eq-tte-1

事件记录的精度取决于采集系统的频率，以 20 MHz 的晶振为例，其时钟间隔 $\delta t = 50\ \mathrm{ns}$，因此式 (eq-tte-1) 可以表示为：

\begin{align} D_\mathrm{TTE}: \{m_n,n=1\dots N\} \end{align}

eq-tte-2

其中 $m_n$ 为时间刻度（time ticks），其代表的时间与总的刻度数为：

\begin{align} t_m &= m\delta t \\ M &= \frac{T}{\delta t} \end{align}

进一步，式 (eq-tte-1) 可以使用可观测量 $X_m$ 表示为：

\begin{align} D_\mathrm{TTE}: X_m = \begin{cases} 0, \mathrm{no\ event\ during\ tick\ m} \\ 1, \mathrm{otherwise} \end{cases} \end{align}

其概率可以表示为：

\begin{equation}\begin{split} P(X_m=0|\Lambda) &= e^{-\Lambda} = p_0 \\ P(X_m=1|\Lambda) = 1 - p_0 &= p_1 = e^{-\Lambda}\sum\limits_{k=1}^{\infty}\frac{\Lambda^k}{k!} \\ \Lambda &= \lambda\delta t \end{split}\end{equation}

由于时间间隔 $\delta t$ 极短，所以当 $\lambda$ 不太大时，在一个时间间隔内多于 1 个入射事例的可能性相当低，则有：

\begin{align} p_1 = 1 - e^{-\lambda\delta t} \approx \lambda\delta t = \Lambda \end{align}

eq-tte-p1

Binned Data¶

Binned Data 是已经分箱数据，将时间区间 $T$ 均分为 $M$ 个区间，$X_m$ 表示第 $m$ 个区间的事例计数，则有：

\begin{align} D_\mathrm{BIN}: \{X_m,m=1,\dots,M\} \\ P(X_m|\Lambda) = \frac{\Lambda^{X_m}}{X_m!}e^{-\Lambda} \end{align}

Time-to-Spill Data¶

Time-to-Spill（TTS）是指为了减少数据量，而每隔 $S$ 个事例记录一个事例数据的模式，可表示为：

\begin{align} D_\mathrm{TTS}: \{\tau_n,n=1,\dots,N-1\} \end{align}

$\tau_n$ 是第 $n$ 个和第 $n+1$ 个记录的事例之间的时间间隔（以刻度数表示，即实际时间间隔为 $\tau_n\delta t$），其分布满足：

\begin{align} P(\tau|\Lambda) = \frac{\Lambda^S}{\Gamma(S)}\tau^{S-1}e^{-\Lambda\tau} \end{align}

混合模式¶

实际观测过程中，同一事件的数据模式可能会发生改变。例如计数率太高时可能从 TTE 模式更改为 BIN 或者 TTS 模式，但只要选择的块评价函数对于不同数据模式而言是一致的，则可以处理任意混合模式的数据。

曝光¶

在观测过程中，受到各种因素的影响，仪器的相应并不是一致的，例如辐射源的强度可能会改变探测器的探测效率，天气和环境的改变也可能会造成遮挡或噪声。这些因素的综合影响被称为曝光（Exposure），可以直接从整个探测系统（包含环境）的性质出发进行计算得到，或者通过仪器的刻度与标定获得。在实际中我们可以将第 $n$ 个数据点的曝光表示为：

\begin{align} 0\le e_n \le 1 \end{align}

在实际计算中，需要将曝光的影响进行归一化，即将计数率乘以一个因子 $\frac{1}{e_n}$。对于 TTE 数据我们得到的只是离散的时间序列，为了将其转化为计数率，需要除以事例间隔长度，具体而言对于第 $n$ 个事例，其计数率 $C_n$ 可以表示为：

\begin{align} C_n &= \frac{1}{e_n\Delta t_n} \\ \Delta t_n &= (t_{n+1}-t_{n-1}) / 2 \end{align}

可以证明 $C_n$ 满足一个仅与第 $n$ 个事例时的真实计数率相关的分布。

近似算法¶

常数泊松分布模型¶

假定在时间区间 $T$ 内，辐射源的强度恒定不变，用单位时间的光子数 $\lambda > 0$ 表示。将观测时段以固定时间间隔 $\Delta T$ 划分为若干子区间，每个子区间上的计数 $k$ 均满足泊松分布，即：

\begin{equation}\begin{split} P(k|\Lambda) &= \frac{\Lambda^ke^{-\Lambda}}{k!} \\ \Lambda &\equiv \lambda\Delta T \end{split}\end{equation}

使用 $\mathscr{M}(\Lambda, T)$ 表示该模型类。除了常数模型外，还有线性模型与指数模型等，其表达式分别为：

\begin{align} \lambda(t) &= \lambda_0[1+a(t-t_0)] \\ \lambda(t) &= \lambda_0 e^{a(t-t_0)} \end{align}

其中 $\lambda_0$ 是基准强度，$t_0$ 是基准时间，$a$ 是强度的变化速率。

TTE Data¶

由于不同入射事例之间是彼此独立的，则 $D_\mathrm{TTE}$ 的联合概率可以表示为每一刻度上概率的乘积：

\begin{equation}\begin{split} P[D_\mathrm{TTE}|\mathscr{M}(\Lambda, T)] &= \prod\limits_{m=1}^MP(X_m|\Lambda) \\ &= p_1^Np_0^{M-N} \end{split}\end{equation}

eq-tte-join

可以发现与参数 $T$ 无关，模型类可简写为 $\mathscr{M}(\Lambda)$。由二项分布可知当 $p_1$ 满足：

\begin{align} \frac{N}{M+1}\le p_1 \le \frac{N+1}{M+1} \end{align}

时，式 (eq-tte-join) 有最大值，由于 $N$ 和 $M$ 都较大，不妨取 $p_1=\frac{N}{M}$。则由式 (eq-tte-p1) 有：

\begin{align} \lambda &= -\frac{1}{\delta t}\ln\left(1-\frac{N}{M}\right) \\ &\approx \frac{1}{\delta t}\frac{N}{M} \end{align}

eq-tte-lambda-max-1

eq-tte-lambda-max-2

使用 $p_1$ 代替 $\Lambda$ 作为模型参数，设参数 $p_1$ 的先验分布为：

\begin{align} P(p_1|\mathscr{M}) = \begin{cases} 1, 0 \le p_1 \le 1 \\ 0, \mathrm{otherwise} \end{cases} \end{align}

代入式 (eq-global-likelihood) 中的全局似然函数，可得：

\begin{equation}\begin{split} P(D_\mathrm{TTE}|\mathscr{M}) &= \int P[D_\mathrm{TTE}|\mathscr{M}(p_1)]P(p_1|\mathscr{M})\mathrm{d}p_1 \\ &= \int_0^1 p_1^N(1-p_1)^{M-N}\mathrm{d}p_1 \\ &= B(N+1, M-N+1) \end{split}\end{equation}

其中 $B(x,y)$ 为贝塔函数，可以用伽马函数表示为：

\begin{align} B(x,y) = \frac{\Gamma(x)\Gamma(y)}{\Gamma(x+y)} \end{align}

因此可得：

\begin{equation}\begin{split} \mathscr{L}(\mathscr{M}|D_\mathrm{TTE}) &= \frac{\Gamma(N+1)\Gamma(M-N+1)}{\Gamma(M+2)} \\ &= \frac{N!(M-N)!}{(M+1)!} \end{split}\end{equation}

eq-tte-likelihood

BIN Data¶

同上可得 $D_\mathrm{BIN}$ 的联合概率为：

\begin{align} P[D_{BIN}|\mathscr{M}(\Lambda)] = \frac{\Lambda^Ne^{-M\Lambda}}{\prod\limits_{m=1}^MX_m!},\ N = \sum\limits_{m=1}^MX_m \end{align}

eq-bin-joint

将式 (eq-bin-joint) 中的分子对 $\Lambda$ 求导可得：

\begin{align} \frac{\mathrm{d}(\Lambda^Ne^{-M\Lambda})}{\mathrm{d}\Lambda} = N\Lambda^{N-1}e^{-M\Lambda}-M\Lambda^Ne^{-M\Lambda} \end{align}

令其等于 0 则可得最大概率的 $\Lambda$ 取值为：

\begin{align} \Lambda = \frac{N}{M} \end{align}

eq-bin-lambda-max

与式 (eq-tte-lambda-max-2) 中的近似值一致。设参数 $\Lambda$ 的先验分布为：

\begin{align} P(\Lambda|\mathscr{M}) = \begin{cases} \frac{e^{-\Lambda}}{1 - e^{-C}}, 0 \le\Lambda\le C \\ 0, \mathrm{otherwise} \end{cases} \end{align}

eq-bin-prior

可以证明该分布等价于 $p_0$ 在 $[e^{-C},1]$ 上的均匀分布。此处的 $C$ 可以近似用 $\frac{T}{Mt_d}$ 表示，其中 $t_d$ 表示仪器的死时间。对于特定的一组数据，不同模型中式 (eq-bin-joint) 中的分母不会发生改变，可以略去。当考虑 $C\to\infty$ 时，则有：

\begin{equation}\begin{split} \mathscr{L}(M|D_\mathrm{BIN}) = \int_0^\infty\Lambda^Ne^{-(M+1)\Lambda}\mathrm{d}\Lambda = \frac{\Gamma(N+1)}{(M+1)^{N+1}} \end{split}\end{equation}

eq-bin-likelihood

TTS Data¶

同上可得 $D_\mathrm{TTS}$ 的联合概率为：

\begin{align} P[D_\mathrm{TTS}|\mathscr{M}(\Lambda)] = \left[\frac{\Lambda^S}{\Gamma(S)}\right]^{N-1}\left(\prod\limits_{n=1}^{N-1}\tau_n\right)^{S-1}e^{-\Lambda M},\ M=\sum\limits_{n=1}^{N-1}\tau_n \end{align}

同理可求得最大概率处 $\Lambda$ 的取值为：

\begin{align} \Lambda = \frac{SN}{M} \end{align}

eq-tts-lambda-max

结合式 (eq-bin-prior)，当 $C\to\infty$ 时可求得似然函数为：

\begin{align} \mathscr{L}(\mathscr{M}|D_\mathrm{TTS}) = \frac{\left(\prod_{n=1}^{N-1}\tau_n\right)^{S-1}}{\Gamma(S)^{N-1}}\frac{\Gamma[S(N-1)+1]}{(M+1)^{S(N-1)+1}} \end{align}

eq-tts-likelihood

总结¶

在常数泊松分布模型下，不同数据情况下的似然度仅取决于时间区间 $T$ 的长度（以刻度数 $M$ 表示）和区间内的事例数目 $N$，因此可以写成：

\begin{align} \mathscr{L}(\mathscr{M}_1|D) = \phi_D(N,M) \end{align}

其中 $\mathscr{M}_1$ 代表常数泊松分布模型。不同 $D$ 情况下的似然函数由式 (eq-tte-likelihood), (eq-bin-likelihood), (eq-tts-likelihood) 给出：

\begin{align} \phi_\mathrm{TTE} &= \frac{\Gamma(N+1)\Gamma(M-N+1)}{\Gamma(M+2)} = \frac{N!(M-N)!}{(M+1)!} \\ \phi_\mathrm{BIN} &= \frac{\Gamma(N+1)}{(M+1)^{N+1}} \\ \phi_\mathrm{TTS} &= \frac{\left(\prod_{n=1}^{N-1}\tau_n\right)^{S-1}}{\Gamma(S)^{N-1}}\frac{\Gamma[S(N-1)+1]}{(M+1)^{S(N-1)+1}} \end{align}

分段常数泊松分布模型¶

在未分段的常数泊松分布模型上改进，将整个时间区间分为若干子区间，每一个子区间内都是恒定的泊松分布模型。考虑最简单的二分段模型，可以使用 $\mathscr{M}_2(\Lambda_1,\Lambda_2,t_\mathrm{cp})$ 表示，该模型可以表示为两个简单模型的复合：

\begin{align} \mathscr{M}_2(\Lambda_1,\Lambda_2,t_\mathrm{cp}) = \mathscr{M}_1(\Lambda_1,M_1)\otimes\mathscr{M}_1(\Lambda_2,M_2) \end{align}

其中 $t_\mathrm{cp}$ 表示两段的分割点（Change Point），在该点处泊松分布的参数从 $\Lambda_1$ 改变为 $\Lambda_2$，前后两段的计数分别为 $N_1$ 与 $N_2$，长度分别为 $M_1$ 与 $M_2$。由于前后两个过程是独立的，因此复合模型的概率可以表示为两个子模型的乘积：

\begin{align} P[D|\mathscr{M}_2(\Lambda_1,\Lambda_2,t_\mathrm{cp})] = P[D_1|\mathscr{M}_1(\Lambda_1,M_1)]\times P[D_2|\mathscr{M}_1(\Lambda_2,M_2)] \end{align}

进而可得到全局似然函数：

\begin{equation}\begin{split} \mathscr{L}(\mathscr{M}_2|D) = &\int\mathrm{d}t_\mathrm{cp}\int\mathrm{d}\Lambda_1\int\mathrm{d}\Lambda_2 P_\mathrm{cp}(t_\mathrm{cp}) \\ &\times P[D_1|\mathscr{M}_1(\Lambda_1,M_1)]P_\Lambda(\Lambda_1|\mathscr{M}_1) \\ &\times P[D_2|\mathscr{M}_1(\Lambda_2,M_2)]P_\Lambda(\Lambda_2|\mathscr{M}_1) \end{split}\end{equation}

积分后的表达式应当与参数 $\Lambda_1$， $\Lambda_2$ 和 $t_\mathrm{cp}$ 无关。两个子模型之间应当是无关的，因此可分离变量为：

\begin{align} \psi(t_\mathrm{cp}) = \phi(N_1,M_1)\phi(N_2,M_2) \\ \mathscr{L}(\mathscr{M}_2|D) = \int P_\mathrm{cp}(t_\mathrm{cp})\psi(t_\mathrm{cp})\mathrm{d}t_\mathrm{cp} \end{align}

实际上各种数据中时间点都是离散的，因此 $\mathrm{d}t_\mathrm{cp}$ 的积分可以使用求和代替。考虑 TTE 情况下的数据，可以使用时间刻度来表示分割点。在没有事例记录的时间点进行分割是没有意义的，因此分割点应当在式 (eq-tte-2) 表示的数据点中，使用 $m_\mathrm{cp} = m_{n_\mathrm{cp}}$ 表示，则有：

\begin{equation}\begin{split} N_1 &= n_\mathrm{cp} \\ N_2 &= N - N_1 \\ M_1 &= m_\mathrm{cp} \\ M_2 &= M - M_1 \end{split}\end{equation}

设 $m_\mathrm{cp}$ 的先验分布为：

\begin{align} P(m_\mathrm{cp}|\mathscr{M}_2) = \begin{cases} \frac{1}{N}, m_\mathrm{cp}\in D_\mathrm{TTE} \\ 0, \mathrm{otherwise} \end{cases} \end{align}

则似然函数可以表示为：

\begin{align} \psi(n_\mathrm{cp}) &= \phi(N_1,M_1)\phi(N_2,M_2)\\ \mathscr{L}(\mathscr{M}_2|D) &= \frac{1}{N}\sum\limits_{n_\mathrm{cp}=1}^N\psi(n_\mathrm{cp})\Delta t_{n_\mathrm{cp}} \end{align}

eq-piecewise-psi

eq-piecewise-likelihood

式 (eq-piecewise-likelihood) 前的常数 $\frac{1}{N}$ 可忽略，而 $\Delta t_{n_\mathrm{cp}}$ 是相邻两个事例间的时间间隔，一般而言是一个可忽略的小修正项。

分块数目¶

根据式 (eq-odds-ratio) 有：

\begin{align} \mathscr{O}_{21} = \frac{J(\mathscr{M}_2|D)}{J(\mathscr{M}_1|D)} = \rho\frac{\mathscr{L}(\mathscr{M}_2|D)}{\mathscr{L}(\mathscr{M}_1|D)} \end{align}

当先验概率比值 $\rho=1$ 时，该比值被称为贝叶斯因子（Bayes Factor，BF）。如果分段模型更占优，则由式 (eq-piecewise-psi) 计算最可能的分段点的位置，$\Lambda$ 的最可几取值可由式 (eq-tte-lambda-max-1)，(eq-tte-lambda-max-2)，(eq-bin-lambda-max)，(eq-tts-lambda-max) 确定。对于具有 $N$ 个分割点的模型，可以逐一求解分割点位置，也可以使用迭代的方法进行近似计算。

迭代方法对整体区间分为两段，如果 BF 更支持分段模型，则对分段后的子区间继续同样的过程，该方法得到的分段方案与直接计算得到的方案十分接近。结束条件可以是所有的子区间 BF 都更支持不分段模型，但是大量的数据显示绝大部分的 $\mathscr{O}_{21}$ 都比 1 大一点，因此会存在区间划分过细的问题。

一个方法是设置子区间允许的最小事例数目以保持良好的统计性，当事例数目达到或小于这一设定值时不再进行分段；另一种方案是改变模型的先验概率值使其更倾向于不分段的模型，即调整 $\mathscr{O}_{21}$ 中 $\rho$ 的大小。

DP 算法¶

块评价函数¶

块评价函数用于评价划分出的块的优劣，以此为基础解决分块的最优化问题。贝叶斯块算法要求块评价函数具有可加性（block-additive），即分组的整体评价是各个块的评价函数值之和。对于某个块评价函数 $f(B_k)$ 和某种分组 $P$，分组的整体评价为：

\begin{align} F(P)=\sum\limits_{k=1}^{N_\mathrm{b}}f(B_k) \end{align}

其中 $B_k$ ​代表第 $k$ 个块，$N_\mathrm{b}$ 为分组 $P$ 中分块的数目。特定的分组方式被记为 $P_\mathrm{cps}$，$P(k)$ 表示前 $k$ 个数据点的分组方式。贝叶斯块算法的目标是选取恰当的分组方式，使得整体评价最优，可以表示为：

\begin{align} \max\limits_{\mathrm{cps}}F(P_\mathrm{cps}) \end{align}

实现步骤¶

假设 $P_\mathrm{opt}(k)$ 为前 $k$ 个点的最佳分组，根据 $P_\mathrm{opt}(k+1)$ 的定义，其最后一个块的右边界一定为第 $k+1$ 个点。定义 $r$ 满足 $P_\mathrm{opt}(k+1)$ 的最后一个块的左边界为第 $r$ 个点，则 $r$ 共有 $k+1$ 种可能，即 $1,2,\dots,k+1$，最后一个块的评价表示为 $f(r)$。

可以证明对于所有以 $r$ 为最后一个块的左边界的 $P(k+1)$ 分组的集合，仅有 $P_\mathrm{opt}(r-1)$ 与这最后一个块组成的方案才可能为最佳分组：

\begin{align} F[P_\mathrm{cps}(r-1)] &\le F[P_\mathrm{opt}(r-1)] \nonumber \\ F[P_\mathrm{cps}(r-1)] + f(r) &\le F[P_\mathrm{opt}(r-1)] + f(r) \end{align}

由于块评价函数具有可加性，因此:

\begin{equation}\begin{split} F[P^r(k+1)] = f(r) + \begin{cases} 0 &; r=1\cr F[P_\mathrm{opt}(r-1)] &; r=2,3,...,R+1 \end{cases} \end{split}\end{equation}

​当 $r=1$ 时，所有数据点都在同一个块中，不存在 $r-1$ 对应的块评价函数。遍历 $k+1$ 种 $r$ 的取值可能，然后取 $r_\mathrm{opt}$ 满足：

\begin{align} r_\mathrm{opt}=\mathop{\mathrm{argmax}}\limits_{r}F[P^r(k+1)] \end{align}

结合 $P_\mathrm{opt}(r_\mathrm{opt}-1)$ 可得 $P_\mathrm{opt}(k+1)$，进而得到 $P_\mathrm{opt}(k), k=1,2,\dots,N$，最优的分组方案即 $P_\mathrm{opt}(N)$。

分块数目¶

$N_\mathrm{b}$ 的先验分布，如：

\begin{align} P(N_\mathrm{b}) = P_0\gamma^{N_\mathrm{b}} \end{align}

其中 $0 \le N_\mathrm{b} \le N$，$0 \le \gamma \le 1$ 是唯一的参数，$P_0$ 和 $\overline{N_\mathrm{b}}$ 可以计算得出：

\begin{align} P_0 &= \frac{1 - \gamma}{1 - \gamma^{N+1}} \\ \overline{N_\mathrm{b}} &= \frac{N\gamma^{N+1}+1}{\gamma^{N+1}-1}+\frac{1}{1-\gamma} \end{align}

定义误报率 $p_1$ 为，正确率 $p_0$ 为：

\begin{align} p_0 \equiv 1 - p_1 \end{align}

未完待续

参考资料¶

1. Scargle J D. Studies in astronomical time series analysis. V. Bayesian blocks, a new method to analyze structure in photon counting data[J]. The Astrophysical Journal, 1998, 504(1): 405.
2. Scargle J D, Norris J P, Jackson B, et al. Studies in astronomical time series analysis. VI. Bayesian Block representations[J]. The Astrophysical Journal, 2013, 764(2): 167.
3. bayesian_blocks