title: Gap Statistic Method
body: <pre><code>    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&gt;0:
              opt_k = i+10
              break
      print diff
      plt.plot(np.linspace(10,19,10,True),gap)
      plt.show()
  </code></pre>
  
  <p>Here I am trying to implement the <strong>Gap Statistic method</strong> for determining the optimal number of clusters. But the problem is that every time I run the code I get a different value for k.
  What is the solution to the problem?
  How can the value of optimal k differ for the same data?</p>
  
  <p><em>I have stored the data in a <code>.mat</code> file beforehand and I am passing it as an argument via terminal</em></p>
  
  <p>I am looking for the smallest value of k for which <code>Gap(k)&gt;= Gap(k+1)-s(k+1)</code> where <code>s(k+1) = sd(k+1)*square_root(1+(1/B))</code> where sd is the standard deviation of the reference distribution and B is the number of copies of Monte Carlo sample</p>
  
  <p>Otherwise stated, I am searching for the value of k for which </p>
  
  <p><code>s(k+1)-Gap(k+1)+Gap(k)&gt;=0</code></p>