Статистика разрыва со стандартной ошибкой 1

Question

Я реализовал Kmeans с помощью команды Scikit Learn и попробовал Elbow и Silhoutte Coefficient, чтобы найти оптимальный K. Я планирую использовать статистику гэпа для дальнейшей проверки моих результатов.

def optimalK(data, nrefs=3, maxClusters=15):


gaps = np.zeros((len(range(1, maxClusters)),))
resultsdf = pd.DataFrame({'clusterCount':[], 'gap':[]})
for gap_index, k in enumerate(range(1, maxClusters)):

    # Holder for reference dispersion results
    refDisps = np.zeros(nrefs)

    for i in range(nrefs):

        # Create new random reference set
        randomReference = np.random.random_sample(size=data.shape)

        # Fit to it
        km = KMeans(k)
        km.fit(randomReference)

        refDisp = km.inertia_
        refDisps[i] = refDisp

    km = KMeans(k)
    km.fit(data)

    origDisp = km.inertia_

    # Calculate gap statistic
    gap = np.log(np.mean(refDisps)) - np.log(origDisp)

    # Assign this loop's gap statistic to gaps
    gaps[gap_index] = gap

    resultsdf = resultsdf.append({'clusterCount':k, 'gap':gap}, ignore_index=True)

return (gaps.argmax() + 1, resultsdf)  

Однако мои графики для статистики разрыва растут, поэтому оптимальное количество кластеров всегда является конечной точкой для моего диапазона кластеров. Предположим, я определяю диапазон кластеров от 1 до 10, тогда оптимальным будет 10.

Согласно веб-сайтам в Интернете и оригинальной статье, обходной путь заключается в реализации стандартной ошибки 1, в которой

GAP(K)> GAP(K+1)- S(K+1)

Может кто-нибудь объяснить мне, как реализовать это в приведенном выше коде? Я не знаю, как рассчитать S(k+1), поскольку это включает в себя поиск стандартного отклонения эталонного распределения.

s(k+1) = sd(k+1)*square_root(1+(1/B))

B - количество копий образцов Монте-Карло. Я смотрю на разные сайты, но кажется, что они не реализовали статистику разрывов со стандартной ошибкой 1.

Answer 1

def gap_stat(data,label):
    k = len(np.unique(label))
    n = data.shape[0]
    p = data.shape[1]
    D_r = []
    C_r = []
    for label_number in range(0,k):
        this_label_index = np.where(label==label_number)[0]
        temp_sum = 0
        pairwise_distance_matrix = 
    euclidean_distances(data[this_label_index],squared=True)
        D_r.append(np.sum(pairwise_distance_matrix))
        C_r.append(float(len(this_label_index)))

    W_r = np.sum(np.asarray(D_r)/(2*np.asarray(C_r)))
    gap_stats = np.log(float(p*n)/12)-(2/float(p))*np.log(k)- 
    np.log(W_r)
    return(gap_stats)