python 检测曲线拐点 的算法 对噪声不敏感 python求曲线拐点

1 最优分割法简介

最有分割法是对有序样品的一种聚类方法。当样品是按顺序排列，在分类中不允许打破样品的顺序。即

，对 个有序样品进行分割，就可能有

种划分方法，这每一种分法称为一种分割。在所有的这些分割中，找到一种分割法，这种分割法使得各段内样品之间的差异最小，而各段之间的差异最大。这种对 个样品分段并使组内离差平方和最小的分割方法，称为最优分割法。段内数值变化可用变差或者极差来表示，比如样品段

变差： 

，

表示样品段内样本间的差异情况，越小表示段内各样品之间差异性越小；因此，把 称为段直径。

2 Fisher最优分割法简介

Fisher最优分割法是用离差平方和来表示同类样本之间的差异程度，通过简便的计算步骤和作图，确定最优分类数，使同类样本间的差异最小，各类别样本间的差异最大，并用F检验法检验最优分类数的合理性。

设某一类 包含的样品有

，记为

，该类的均值向量 为  

， ，

表示这一类的直径，常用的直径有

，当是单变量时，也可定义直径为­

， 是中位数。

3 Fisher最优分割递推公式

用

表示将n 个有序样品分为 类的一种分法，定义损失函数为

，当 和 固定时，

越小表示各类的离差平方和越小，因此要寻找一种分法

，使分类损失函数

达到最小，记为

。递推公式：

             

公式含义：如果要找到n 个样品分为 个类的最优分割，应建立在将

个样品分为k-1 类的最优分割的基础上。

4 Fisher最优分割计算流程

若分类数 已知，求分类法

使它在损失函数意义下达最小，其求法如下：首先找分点 使递推公式达极小，即

。于是得第k类

。然后找 ，使它满足

得到第

类。类似的方法一次可以得到所有类

,这就是所求的最优解，即：

。总之，为了求最优解，主要是计算

和

}

5 Fisher最优分割-Python完整代码：

import numpy as np
def G(X, i, j):
    #tex:
    #$$ \overline{x_G}=\frac{1}{j-i+1}\sum^{j}_{t=i}X_t $$
    xg = 0
    for t in range(i-1, j):
        xg += X[t]
    xg /= (j - i + 1)
    return xg
def D(X, i, j):
    #tex:
    #$$ D(i,j)=\sum^{j}_{t=i}(X_t-\overline{x_G})'(X_t-\overline{x_G}) $$
    d = 0
    xg = G(X, i, j)
    for t in range(i-1, j):
        d += (X[t] - xg) ** 2
    return d
def Lq(X, N, K):
    #tex:
    #$$ L(b(n,k))=\sum^k_{i=1}D(i_t,i_{t+1}-1) $$
    #$$ \begin{equation} \left\{ \begin{array}{lr} L(p(n,2))=1 \\ L(p(n,k))=\underset{k \le j \le n}{min} D(1,j-1)+D(n,j) \end{array} \right. \end{equation} $$
    l = np.zeros((N+1, K+1))
    for n in range(1, N+1):
        l[n, 1] = D(X, 1, n)
    for k in range(2, K+1):
        for n in range(k, N+1):
            l[n, k] = min([l[j-1, k-1] + D(X, j, n) for j in range(k, n+1)])
    return l
def B(Lp):
    N = Lp.shape[0] - 1
    K = Lp.shape[1] - 1
    Bl = np.zeros(K+1)
    for k in range(1, K):
        Bl[k] = Lp[N, k] / Lp[N, k+1]
    return Bl
def split_class(X, K, l=None):
    if l is None:
        l = Lq(X, len(X), K * 2 + 1)
        Bl = B(l)
        for i in range(1, K):
            if Bl[i]-1 <= 0.05:
                K = i+1
                l = Lq(X, len(X), K * 2 + 1)
                break
    N = len(X)
    if K == 1:
        return []
    else:
        dist = np.array([abs(l[t-1, K-1] + D(X, t, N) - l[N, K]) for t in range(K, N+1)])
        t = list(range(K, N+1))[int(np.argmin(dist))]
        r = split_class(X[:t-1], K-1, l)
        r.append(t - 1)
        return r

直接调用split_class函数即可

6 Fisher最优分割-确定分类数K

分类个数的确定如果能从实际问题中首先确定当然最好，如果不能则可以从

随k的变化趋势图中找到拐点处，作为确定k的依据。当曲线拐点很平缓时，可选择的k较多，这时需要用其他方法来确定，如均方比法和特征根法。稳定分类数可以这样定义: 当分类数小于此值时, 

较大, 类内样本间的差别仍然显著, 需要继续再细分, 否则分类选取时, 代表性和信息量不能保证; 当分类数大于此值时, 类内样本间的差别已经不显著, 再细分的意义不大。稳定分类数代表了

值由较大值向较小值变化的转折点。按分类选取的原则在每一类中选取一个代表点, 则在优化布点时, 分类数应不小于稳定分类数, 才能保证优选点的代表性。