转载请标明出处：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; heta) = f(x_1, ..., x_n; heta)$：1）当$X$是连续型的，则$f(x_1,...,x_n; heta)$是其联合概率密度，如果$X$是简单随机样本，即$X_1,...,X_n$是独立同分布的，概率函数可以表示成连乘形式$f(x_1,...,x_n; heta) = f(x_1; heta) cdot ... cdot f(x_n; heta) $；2）当$X$是离散型的，则有$f(x_1,...,x_n; heta) = P(X_1 = x_1, ..., X_n = x_n ; ; heta) $，同样地，如果$X$是简单随机样本，有连乘的形式$f(x_1,...,x_n; heta) = P(X_1 = x_1; heta) cdot ...cdot P(X_n = x_n ; heta) $.
定义 1：设样本$X=(X_1, X_2, ..., X_n )$有概率函数$f(x; heta) = f(x_1, ..., x_n; heta) $，其中参数$ heta in Theta $. 当视$x$为常量，$ heta$为定义在参数空间$Theta$上的自变量时，称函数$f(x, heta)$为似然函数。
例 1：设有样本$X=(X_1, X_2, ..., X_n )$，$X_1,...,X_n sim N(mu, sigma^2)$相互独立。此时参数为$ heta = (mu, sigma^2)$，参数空间$Theta = { (mu, sigma^2) in mathbb R^2 : sigma^2 ge 0 } $，根据定义 1上方的分析可计算似然函数：$$ egin{equation} egin{split} f(x; heta) & = 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}{2sigma^2} ight ) end{split} end{equation} $$从上述表达式易知，当样本观测值$(x_1, ..., x_n)$给定时，似然函数$f(x; heta)$因取不同参数$ heta = (mu, sigma^2)$时而出现不同的结果。
然而如果$X_1,...,X_n sim N(mu, sigma^2) $但不独立，此时因此样本变量之间不独立，参数$ heta $中将除$mu $和$sigma^2 $以外还包含两两样本变量之间的协方差，假设协方差矩阵为$Sigma $，似然函数如下式：$$ egin{equation} f(x; heta) = frac{1}{(2pi)^{n/2} |Sigma|^{1/2}}exp left ( - frac12 (vec x - mu)^T Sigma^{-1} (vec x - mu) ight ) end{equation} $$其中向量$vec x = (x_1,...,x_n)^T $.
从定义可知，样本分布的函数形式是已知的，只是因为参数$ heta $还是未知的，一旦确定参数$ heta $的值，即可完全给定样本分布。反过来看，现在已经有了一组结果（样本观测值）$(x_1, ..., x_n) $，任意的参数$ heta $值都会反馈到似然函数的输出，该输出正比于导致这一组结果的可能性。如果有参数$ heta^* $代入似然函数计算后得到最大输出，则有理由相信已知的一组样本观测值最有可能服从参数为$ heta^* $的样本分布，这就是极大似然估计的思想：
定义 2：若$hat heta(X) $是一个统计量，满足条件$$ egin{equation} f(x;hat heta(x)) = sup_{ heta in Theta} f(x; heta) quad (x in mathfrak X) end{equation} $$其中，$mathfrak X $是样本空间，则称$hat heta(X) $是$ heta $的极大似然估计。若待估函数是$g( heta) $，则称$g(hat heta(X)) $为极大似然估计。
在统计上，把凡是由样本算出的量称为统计量，最为常用的统计量是样本均值$$ egin{equation} ar X = frac1n sum_{i=1}^n X_i end{equation} $$和样本方差$$ egin{equation} S^2 = frac1{n-1} sum_{i=1}^n (X_i - ar X)^2 end{equation} $$样本$X=(X_1, X_2, ..., X_n ) $可能取值的全体的集合称为样本空间。定义中参数的函数$g( heta) $一般指参数$ heta $中的某个或某几个分量，如在例 1中，参数$ heta = (mu, sigma) $，可以令$g( heta) = mu $，表示只估计参数$ heta $的第1个分量$mu $．
对数似然函数为$ln f(x; heta) $．如果$X=(X_1, X_2, ..., X_n ) $为简单随机样本，似然函数可用连乘形式计算：$$ egin{equation} f(x; heta) = prod_{i=1}^n f(x_i; heta) end{equation} $$此时，对数似然函数有更简单的形式：$$ egin{equation} ln f(x; heta) = sum_{i=1}^n ln f(x_i; heta) end{equation} $$因为对数函数是严格单调增的，(3)式中的极大似然估计$hat heta (X) $可以等价地定义为：$$ egin{equation} ln f(x;hat heta(x)) = sup_{ heta in Theta} ln f(x; heta) end{equation} $$
式(3)和(8)中的极值一般选择下述两个似然方程中的一个进行求解：$$ egin{equation} frac{mathrm d f(x; heta)}{mathrm d heta} = 0 end{equation} $$$$ egin{equation} frac{mathrm d ln f(x; heta)}{mathrm d heta} = 0 end{equation} $$
例 2：设简单随机样本$X=(X_1,...,X_n)sim N(mu, sigma^2) $，参数$ heta = (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 $求偏导得似然方程$$ egin{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 = ar x $$$$sigma^2 = frac 1n sum_{i=1}^n (x_i - ar x)^2$$因此参数的极大似然估计为$hat mu(X) = ar X $，$hat sigma^2 = 1/n sum_{i=1}^n (X_i - ar 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; heta) = f(x_1,...,x_n; p_1,...,p_k) = p_1^{n_1} cdots p_k^{n_k} $$其中参数$ heta = (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; heta) = sum_{i=1}^{k-1} n_i ln p_i + n_k ln (1-sum_{i=1}^{k-1} p_i)$$对参数求偏导得似然方程$$ egin{split} & frac{partial}{partial p_i} ln f(x; heta) = frac{n_i}{p_i} - frac{n_k}{1-sum_{j=1}^{k-1} p_j} = 0, quad iin { 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 = ( ilde x_{i1}, ..., ilde x_{iD_1}, ddot x_{i1},..., ddot x_{iD_2}) $，其中$ ilde 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 = ( ilde X_1,..., ilde X_{D_1}, ddot X_1,...,ddot X_{D_2}) $，这样$mathfrak D $中每个数据对象$vec x_i $都是样本$vec X $的一个观测值。假设随机向量$vec X $中随机变量之间均相互独立，对应于连续型属性的随机变量服从正态分布，即$ ilde X_k sim N(mu_{jk}, sigma^2_{jk}) ; (k=1,...,D_1) $；对应于分类型属性的随机变量服从多项分布，即$ ddot X_k sim M(1;p_{jk1},..., p_{jkepsilon_k} ) ; (k=1,...,D_2) $，其中$p_{jkl} ; (l=1,..., epsilon_k) $表示第$l $种取值在试验中出现的概率。据此可得样本$vec X $的似然函数：$$ egin{equation} f_j(vec x; heta_j) = prod_{k=1}^{D_1} f( ilde x_k; mu_{jk}, sigma^2_{jk}) prod_{k=1}^{D_2} f(ddot x_k; p_{jk1},...,p_{jkepsilon_k}) end{equation} $$得对数似然函数：$$ egin{equation} ln f_j(vec x; heta_j) = sum_{k=1}^{D_1} ln f( ilde x_k; mu_{jk}, sigma^2_{jk}) + sum_{k=1}^{D_2} ln f(ddot x_k; p_{jk1},...,p_{jkepsilon_k}) end{equation} $$根据例 2有$$ egin{equation} ln f( ilde 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}( ilde x_{ik}-mu_{jk})^2 end{equation} $$参数的极大似然估计为$hat mu_{jk} = 1/N_j sum_{i=1}^{N_j} ilde X_{ik} = ar { ilde X}_k $，$hat sigma^2_{jk} = 1/{N_j} sum_{i=1}^{N_j} ( ilde X_{ik} - ar { ilde X}_k)^2 $．根据例 3有$$ egin{equation} ln f(ddot x_k; p_{jk1},...,p_{jkepsilon_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 $的对数似然值：$$ egin{equation} egin{split} L_{C_j} & = ilde 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} ight ) - 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）：$$ egin{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 angle} $，此时新的簇划分对数似然值为$$ egin{equation} L_n = sum_{j in {1,...,J} setminus {s,t}} L_{C_j} + L_{C_{langle s,t angle}} end{equation} $$因此$$ egin{equation} L_c - L_n = L_{C_s} + L_{C_t} - L_{C_{langle s,t angle}} end{equation} $$消除重复项后得$$ egin{equation} L_c - L_n = xi_s + xi_t - xi_{langle s,t angle} end{equation} $$其中$$ egin{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} ight ) end{equation} $$但存在可能引起不可计算的退化情形：1）簇$C_j $中只有一个数据对象；或者2）簇$C_j $的数据对象某些属性下具体完全相同的属性值。这两种情况都存在$ln0 $的无意义值，处理方式为对于式(20)的第1项加个小的扰动，对于第2项设定$0ln0 =0 $，有$$ egin{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} ight ) end{equation} $$在SPSS Modeler的算法说明文档中设定$Delta_k = hat sigma^2_k $，即考虑所有数据时在第$k $个连续型属性下方差的极大似然估计。如下定义簇$C_s $和$C_t $之间的对数似然距离：$$ egin{equation} d(C_s, C_t) = zeta_s + zeta_t - zeta_{langle s,t angle} end{equation} $$
[1] 极大似然估计, 百度百科, http://baike.baidu.com/link?url= YnwADRuBU7AxHoGK_5IMO8Geaf1e5MjpvyUN9CEciwjBPyaa3UbyRdeqaPqa2vGxVD5ZaeKQu_-4sTyYHLBHwjkDJp02tJVLRoclC-E0BCwO4tCT_G2EUcdTxT7a7rN16FTN3LHuNf3BWKPsMCnK3q.
[2] 陈希孺，倪国熙．数理统计学教程．合肥：中国科学技术大学出版社，2009.
[3] Chiu T, Fang D P, Chen J, et al. A robust and scalable clustering algorithm for mixed type attributes in large database environment[C]// ACM SIGKDD International Conference on Knowledge Discovery and Data Mining. 2001:263-268.