转载请标明出处：http://www.cnblogs.com/tiaozistudy/p/log-likelihood_distance.html

本文是“挑子”在学习对数似然距离过程中的笔记摘录，文中不乏一些个人理解，不当之处望多加指正。

对数似然距离是基于统计理论的一种计算簇与簇相异度的方法，最早用于BIRCH层次聚类算法的改进。本文旨在详细介绍对数似然距离的统计学基础、方法思想和计算过程，希望帮助更多地人欣赏它、熟悉它、使用它。

1、极大似然估计（Maximum Likelihood Estimate）

极大似然估计方法是求点估计的另一种方法，1821年首先由德国数学家C. F. Gauss（高斯）提出，但是这个方法通常被归功于英国的统计学家R. A. Fisher（罗纳德·费希尔），他在1922年的论文On the mathematical foundations of theoretical statistics, reprinted in Contributions to Mathematical Statistics (by R. A. Fisher), 1950, J. Wiley & Sons, New York 中再次提出了这个思想，并且首先探讨了这种方法的一些性质，极大似然估计这一名称也是费希尔给的。这是一种目前仍然得到广泛应用的方法^[1]。本节以下内容主要参考《数理统计学教程》^[2]编撰，如果对极大似然估计十分熟悉，可直接跳过本节。

假设大小为$n$的样本$X=(X_1, X_2, ..., X_n )$有概率函数$f(x;\theta) = f(x_1, ..., x_n; \theta)$：1）当$X$是连续型的，则$f(x_1,...,x_n;\theta)$是其联合概率密度，如果$X$是简单随机样本，即$X_1,...,X_n$是独立同分布的，概率函数可以表示成连乘形式$f(x_1,...,x_n;\theta) = f(x_1;\theta) \cdot ... \cdot f(x_n;\theta) $；2）当$X$是离散型的，则有$f(x_1,...,x_n;\theta) = P(X_1 = x_1, ..., X_n = x_n ; \; \theta) $，同样地，如果$X$是简单随机样本，有连乘的形式$f(x_1,...,x_n;\theta) = P(X_1 = x_1; \theta) \cdot ...\cdot P(X_n = x_n ; \theta) $.

定义 1：设样本$X=(X_1, X_2, ..., X_n )$有概率函数$f(x;\theta) = f(x_1, ..., x_n; \theta) $，其中参数$ \theta \in \Theta $. 当视$x$为常量，$\theta$为定义在参数空间$\Theta$上的自变量时，称函数$f(x, \theta)$为似然函数。

例 1：设有样本$X=(X_1, X_2, ..., X_n )$，$X_1,...,X_n \sim N(\mu, \sigma^2)$相互独立。此时参数为$\theta = (\mu, \sigma^2)$，参数空间$\Theta = \{ (\mu, \sigma^2) \in \mathbb R^2 : \sigma^2 \ge 0 \} $，根据定义 1上方的分析可计算似然函数：$$ \begin{equation} \begin{split} f(x;\theta) & = \prod_{i=1}^n f(x_i; \mu, \sigma^2) \\ & = \prod_{i=1}^n \frac{1}{\sqrt{2 \pi}\sigma} \exp \left ( -\frac{(x_i-\mu)^2}{2\sigma^2} \right ) \end{split} \end{equation} $$从上述表达式易知，当样本观测值$(x_1, ..., x_n)$给定时，似然函数$f(x;\theta)$因取不同参数$\theta = (\mu, \sigma^2)$时而出现不同的结果。

然而如果$X_1,...,X_n \sim N(\mu, \sigma^2) $但不独立，此时因此样本变量之间不独立，参数$\theta $中将除$\mu $和$\sigma^2 $以外还包含两两样本变量之间的协方差，假设协方差矩阵为$\Sigma $，似然函数如下式：$$ \begin{equation} f(x;\theta) = \frac{1}{(2\pi)^{n/2} |\Sigma|^{1/2}}\exp \left ( - \frac12 (\vec x - \mu)^T \Sigma^{-1} (\vec x - \mu) \right ) \end{equation} $$其中向量$\vec x = (x_1,...,x_n)^T $.

从定义可知，样本分布的函数形式是已知的，只是因为参数$\theta $还是未知的，一旦确定参数$\theta $的值，即可完全给定样本分布。反过来看，现在已经有了一组结果（样本观测值）$(x_1, ..., x_n) $，任意的参数$\theta $值都会反馈到似然函数的输出，该输出正比于导致这一组结果的可能性。如果有参数$\theta^* $代入似然函数计算后得到最大输出，则有理由相信已知的一组样本观测值最有可能服从参数为$\theta^* $的样本分布，这就是极大似然估计的思想：

定义 2：若$\hat \theta(X) $是一个统计量，满足条件$$ \begin{equation} f(x;\hat \theta(x)) = \sup_{\theta \in \Theta} f(x; \theta) \quad (x \in \mathfrak X) \end{equation} $$其中，$\mathfrak X $是样本空间，则称$\hat \theta(X) $是$\theta $的极大似然估计。若待估函数是$g(\theta) $，则称$g(\hat \theta(X)) $为极大似然估计。

在统计上，把凡是由样本算出的量称为统计量，最为常用的统计量是样本均值$$ \begin{equation} \bar X = \frac1n \sum_{i=1}^n X_i \end{equation} $$和样本方差$$ \begin{equation} S^2 = \frac1{n-1} \sum_{i=1}^n (X_i - \bar X)^2 \end{equation} $$样本$X=(X_1, X_2, ..., X_n ) $可能取值的全体的集合称为样本空间。定义中参数的函数$g(\theta) $一般指参数$\theta $中的某个或某几个分量，如在例 1中，参数$\theta = (\mu, \sigma) $，可以令$g(\theta) = \mu $，表示只估计参数$\theta $的第1个分量$\mu $．

对数似然函数为$\ln f(x;\theta) $．如果$X=(X_1, X_2, ..., X_n ) $为简单随机样本，似然函数可用连乘形式计算：$$ \begin{equation} f(x;\theta) = \prod_{i=1}^n f(x_i;\theta) \end{equation} $$此时，对数似然函数有更简单的形式：$$ \begin{equation} \ln f(x;\theta) = \sum_{i=1}^n \ln f(x_i;\theta) \end{equation} $$因为对数函数是严格单调增的，(3)式中的极大似然估计$\hat \theta (X) $可以等价地定义为：$$ \begin{equation} \ln f(x;\hat \theta(x)) = \sup_{\theta \in \Theta} \ln f(x; \theta) \end{equation} $$

式(3)和(8)中的极值一般选择下述两个似然方程中的一个进行求解：$$ \begin{equation} \frac{\mathrm d f(x;\theta)}{\mathrm d \theta} = 0 \end{equation} $$$$ \begin{equation} \frac{\mathrm d \ln f(x;\theta)}{\mathrm d \theta} = 0 \end{equation} $$

例 2：设简单随机样本$X=(X_1,...,X_n)\sim N(\mu, \sigma^2) $，参数$\theta = (\mu, \sigma^2) $，$x=(x_1,x_2,...,x_n) $是一组来自样本$X $的样本观测值，极大似然估计$\hat \mu(x) $和$\hat \sigma^2(x) $如下计算：

从例 1知，样本的似然函数为(1)，因此$$ \ln f(x; \mu, \sigma^2) = -\frac n2 \ln (2 \pi) - \frac n2 \ln \sigma^2 - \frac 1{2 \sigma^2} \sum_{i=1}^n(x_i-\mu)^2 $$分别对参数$\mu $和$\sigma^2 $求偏导得似然方程$$ \begin{cases} \dfrac{\partial}{\partial \mu} \ln f(x; \mu, \sigma^2) = \dfrac 1{\sigma^2} (\sum \limits _{i=1}^n x_i - n \mu) = 0 \\ \dfrac{\partial}{\partial \sigma^2} \ln f(x; \mu, \sigma^2) = -\dfrac n{2 \sigma^2} + \dfrac 1{2(\sigma^2)^2} \sum \limits _{i=1}^n (x_i-\mu)^2 = 0 \end{cases} $$求解得$$\mu = \frac 1n \sum_{i=1}^n x_i = \bar x $$$$\sigma^2 = \frac 1n \sum_{i=1}^n (x_i - \bar x)^2$$因此参数的极大似然估计为$\hat \mu(X) = \bar X $，$\hat \sigma^2 = 1/n \sum_{i=1}^n (X_i - \bar X)^2 $．

例 3：设有$k $个事件$A_1,...,A_k $两两互斥，其概率$p_1,...,p_k $之和为1，将试验独立地重复$n $次得样本$X = (X_1,...,X_n) $．此处样本$X_i $服从多项分布$M(1; p_1,...,p_k) $．易知$X $的似然函数为：$$ f(x; \theta) = f(x_1,...,x_n; p_1,...,p_k) = p_1^{n_1} \cdots p_k^{n_k} $$其中参数$\theta = (p_1,...,p_k) $，$n_j \; (j=1,...,k) $表示这$n $次试验中事件$A_j $出现的次数。取对数并代入$p_k = 1 - \sum_{i=1}^{k-1} p_i $得$$\ln f(x; \theta) = \sum_{i=1}^{k-1} n_i \ln p_i + n_k \ln (1-\sum_{i=1}^{k-1} p_i)$$对参数求偏导得似然方程$$ \begin{split} & \frac{\partial}{\partial p_i} \ln f(x;\theta) = \frac{n_i}{p_i} - \frac{n_k}{1-\sum_{j=1}^{k-1} p_j} = 0, \quad i\in \{ 1,...,k-1\} \\ \Rightarrow & n_i p_k = n_k p_i, \quad i \in \{ 1,...,k-1\} \\ \Rightarrow & p_k \sum_{i=1}^{k-1} n_i = n_k \sum_{i=1}^{k-1} p_i \\ \Rightarrow & p_k \sum_{i=1}^k n_i = n_k \sum_{i=1}^{k-1} p_i + n_kp_k =  n_k \\ \Rightarrow & p_k = \frac{n_k}{\sum_{i=1}^k n_i} = \frac{n_k}n \end{split} $$将$ p_k $的结果代入$n_ip_k = n_k p_i $，可计算出极值点下所有参数的极大似然估计：$$ \hat p_i(X) = \frac{n_i}n, \quad i=1,...,k $$

2、簇内对数似然函数

本节和下节内容主要参考文献[3]．设数据集$\mathfrak D $中有$N $个数据对象$\{ \vec x_i: i=1,...,N\} $，每个数据对象由$D $个属性刻画，其中有$D_1 $个连续型属性（continuous attribute）和$D_2 $个分类型属性（categorical attribute），不失一般性，可设$\vec x_i = (\tilde x_{i1}, ..., \tilde x_{iD_1}, \ddot x_{i1},..., \ddot x_{iD_2}) $，其中$\tilde x_{ij} $表示第$i $个数据对象在第$j $连续型属性下的属性值，$\ddot x_{ik} $表示第$i $个数据对象在第$k $分类型属性下的属性值，已知第$k $个分类型属性有$\epsilon_k $种可能取值。假设数据集$\mathfrak D $划分成$J $个簇$C_j \; (j=1,...,J) $，其中簇$C_j $中包含$N_j $个数据对象。

将簇$C_j $中的$N_j $个数据对象视为某个$D $维随机向量$\vec X $的样本，$\vec X = (\tilde X_1,...,\tilde X_{D_1}, \ddot X_1,...,\ddot X_{D_2}) $，这样$\mathfrak D $中每个数据对象$\vec x_i $都是样本$\vec X $的一个观测值。假设随机向量$\vec X $中随机变量之间均相互独立，对应于连续型属性的随机变量服从正态分布，即$\tilde X_k \sim N(\mu_{jk}, \sigma^2_{jk}) \; (k=1,...,D_1) $；对应于分类型属性的随机变量服从多项分布，即$ \ddot X_k \sim M(1;p_{jk1},..., p_{jk\epsilon_k} ) \; (k=1,...,D_2) $，其中$p_{jkl} \; (l=1,..., \epsilon_k) $表示第$l $种取值在试验中出现的概率。据此可得样本$\vec X $的似然函数：$$ \begin{equation} f_j(\vec x; \theta_j) = \prod_{k=1}^{D_1} f(\tilde x_k; \mu_{jk}, \sigma^2_{jk}) \prod_{k=1}^{D_2} f(\ddot x_k; p_{jk1},...,p_{jk\epsilon_k}) \end{equation} $$得对数似然函数：$$ \begin{equation} \ln f_j(\vec x; \theta_j) = \sum_{k=1}^{D_1} \ln f(\tilde x_k; \mu_{jk}, \sigma^2_{jk}) + \sum_{k=1}^{D_2} \ln f(\ddot x_k; p_{jk1},...,p_{jk\epsilon_k}) \end{equation} $$根据例 2有$$ \begin{equation} \ln f(\tilde x_k; \mu_{jk}, \sigma^2_{jk}) = -\frac {N_j}2 \ln (2 \pi) - \frac {N_j}2 \ln \sigma^2_{jk} - \frac 1{2 \sigma^2_{jk}} \sum_{i=1}^{N_j}(\tilde x_{ik}-\mu_{jk})^2 \end{equation} $$参数的极大似然估计为$\hat \mu_{jk} = 1/N_j \sum_{i=1}^{N_j} \tilde X_{ik} = \bar {\tilde X}_k $，$\hat \sigma^2_{jk} = 1/{N_j} \sum_{i=1}^{N_j} (\tilde X_{ik} - \bar {\tilde X}_k)^2 $．根据例 3有$$ \begin{equation} \ln f(\ddot x_k; p_{jk1},...,p_{jk\epsilon_k}) = \sum_{l=1}^{\epsilon_k} N_{jkl} \ln p_{jkl} \end{equation} $$其中，$N_{jkl} \; (l=1,...,\epsilon_k) $表示簇$C_j $中在第$k $个分类型属性下第$l $种取值的数据对象个数，极大似然估计为$\hat p_{jkl} = N_{jkl}/N_j \; (l=1,...,\epsilon_k) $．将极大似然估计代入(12)式的对数似然函数中可得簇$C_j $的对数似然值：$$ \begin{equation} \begin{split} L_{C_j} & = \tilde L_{C_j} + \ddot L_{C_j} \\ & = -\frac{N_j}2 \left ( D_1 \Big [\ln(2 \pi) + 1 \Big ] + \sum_{k=1}^{D_1} \ln \hat \sigma^2_{jk} \right ) - N_j \sum_{k=1}^{D_2} \hat E_{jk} \end{split} \end{equation} $$其中$\hat E_{jk} = -\sum_{l=1}^{\epsilon_k} N_{jkl}/N_j \ln (N_{jkl}/N_j) $表示根据估计出的概率计算出的簇$C_j $中第$k $上分类型属性下的信息熵。

如假设数据集中数据在簇内有相同分布，但簇与簇之间有不同分布参数，则可定义簇划分对数似然值（logarithm of the classification likelihood）：$$ \begin{equation} L = \sum_{j=1}^J L_{C_J} \end{equation} $$

3、对数似然距离

假设当前数据对象被划分成$J $个簇：$C_1,...,C_J $，根据式(16)可计算当前簇划分对数似然值$L_c $，如果将$C_s $和$C_t \; (s,t \in \{1,...,J\}) $两个簇合并成一个簇，记为$C_{\langle s,t \rangle} $，此时新的簇划分对数似然值为$$ \begin{equation} L_n = \sum_{j \in \{1,...,J\} \setminus \{s,t\}} L_{C_j} + L_{C_{\langle s,t \rangle}} \end{equation} $$因此$$ \begin{equation} L_c - L_n = L_{C_s} + L_{C_t} - L_{C_{\langle s,t \rangle}} \end{equation} $$消除重复项后得$$ \begin{equation} L_c - L_n = \xi_s + \xi_t - \xi_{\langle s,t \rangle} \end{equation} $$其中$$ \begin{equation} \xi_j = -N_j \left ( \frac12 \sum_{k=1}^{D_1} \ln \hat \sigma^2_{jk} + \sum_{k=1}^{D_2} \hat E_{jk} \right ) \end{equation} $$但存在可能引起不可计算的退化情形：1）簇$C_j $中只有一个数据对象；或者2）簇$C_j $的数据对象某些属性下具体完全相同的属性值。这两种情况都存在$\ln0 $的无意义值，处理方式为对于式(20)的第1项加个小的扰动，对于第2项设定$0\ln0 =0 $，有$$ \begin{equation} \zeta_j = -N_j \left ( \frac12 \sum_{k=1}^{D_1} \ln (\hat \sigma^2_{jk} + \Delta_k) + \sum_{k=1}^{D_2} \hat E_{jk} \right ) \end{equation} $$在SPSS Modeler的算法说明文档中设定$\Delta_k = \hat \sigma^2_k $，即考虑所有数据时在第$k $个连续型属性下方差的极大似然估计。如下定义簇$C_s $和$C_t $之间的对数似然距离：$$ \begin{equation} d(C_s, C_t) = \zeta_s + \zeta_t - \zeta_{\langle s,t \rangle} \end{equation} $$