Skip to main content
 首页 » 编程设计

python之差距统计法

2025年05月04日8yxwkf

    import sys 
    import numpy as np 
    import scipy.io as sio 
    import matplotlib.pyplot as plt 
    from sklearn.cluster import KMeans 
    from sklearn.svm import SVC 
 
    filename = sys.argv[1] 
    datafile = sio.loadmat(filename) 
    data = datafile['bow'] 
    sizedata=[len(data), len(data[0])] 
    gap=[] 
    SD=[]     
    for knum in xrange(10,20): 
            print knum 
            #Clustering original Data 
            kmeanspp = KMeans(n_clusters=knum,init = 'k-means++',max_iter = 100,n_jobs = 1) 
            kmeanspp.fit(data) 
            dispersion = kmeanspp.inertia_ 
            #Clustering Reference Data 
            nrefs = 10 
            refDisp = np.zeros(nrefs) 
            for nref in xrange(nrefs): 
                refdata = np.random.random_sample((sizedata[0],sizedata[1])) 
                refkmeans = KMeans(n_clusters=knum,init='k-means++',max_iter=100,n_jobs=1) 
                refkmeans.fit(refdata) 
                refdisp = refkmeans.inertia_ 
                refDisp[nref]=np.log(refdisp) 
            mean_log_refdisp = np.mean(refDisp) 
            gap.append(mean_log_refdisp-np.log(dispersion)) 
            #Calculating standard deviaiton 
            sd = (sum([(r-m)**2 for r,m in zip(refDisp,[mean_log_refdisp]*nrefs)])/nrefs)**0.5 
            SD.append(sd) 
    SD = [sd*((1+(1/nrefs))**0.5) for sd in SD] 
    #determining optimal k 
    opt_k = None 
    diff = [] 
    for i in xrange(len(gap)-1): 
        diff = (SD[i+1]-(gap[i+1]-gap[i])) 
        if diff>0: 
            opt_k = i+10 
            break 
    print diff 
    plt.plot(np.linspace(10,19,10,True),gap) 
    plt.show() 

我在这里尝试实现差距统计方法来确定最佳簇数。但问题是,每次我运行代码时,我都会得到不同的 k 值。 解决问题的方法是什么? 对于相同的数据,最佳 k 的值有何不同?

我预先将数据存储在 .mat 文件中,并通过终端将其作为参数传递

我正在寻找最小的 k 值,其中 Gap(k)>= Gap(k+1)-s(k+1) 其中 s(k+1) = sd(k+1)*square_root(1+(1/B)) 其中sd是引用分布的标准差,B是蒙特卡洛样本的拷贝数

另外说明,我正在寻找 k 的值,为此

s(k+1)-Gap(k+1)+Gap(k)>=0

请您参考如下方法:

你的模拟有几个问题:

1- sd = (sum([(r-m)**2 for r,m in zip(refDisp,[mean_log_refdisp]*nrefs)])/nrefs)**0.5 

为什么要将 zip 的第二个组件乘以根据原始论文不需要的 nrefs。

2-

 if diff>0: 
        opt_k = i+10 
        break 

如果 diff>0 你想要 diff>=0 因为平等可以发生 关于为什么每次都会得到不同数量的聚类,正如人们所说,这是蒙特卡洛模拟,因此可能存在随机性,而且它还取决于您要聚类的内容和数据集。我建议您针对 Silhouette 和 Elbow 测试您的算法,以更好地了解集群的数量。