code: public static <T> GaussianRadialBasis[] learnGaussianRadialBasis(T[] x, T[] centers, Metric<T> distance, int p) {
          if (p < 1) {
              throw new IllegalArgumentException("Invalid number of nearest neighbors: " + p);
          }
          
          int k = centers.length;
          CLARANS<T> clarans = new CLARANS<>(x, distance, k, Math.min(100, (int) Math.round(0.01 * k * (x.length - k))));
          System.arraycopy(clarans.medoids(), 0, centers, 0, k);
  
          p = Math.min(p, k-1);
          double[] r = new double[k];
          GaussianRadialBasis[] rbf = new GaussianRadialBasis[k];
          for (int i = 0; i < k; i++) {
              for (int j = 0; j < k; j++) {
                  r[j] = distance.d(centers[i], centers[j]);
              }
              
              Arrays.sort(r);
              double r0 = 0.0;
              for (int j = 1; j <= p; j++) {
                  r0 += r[j];
              }
              r0 /= p;
              rbf[i] = new GaussianRadialBasis(r0);
          }
  
          return rbf;
      }
docs: Learns Gaussian RBF function and centers from data. The centers are
  chosen as the medoids of CLARANS. The standard deviation (i.e. width)
  of Gaussian radial basis function is estimated by the p-nearest neighbors
  (among centers, not all samples) heuristic. A suggested value for
  p is 2.
  @param x the training dataset.
  @param centers an array to store centers on output. Its length is used as k of CLARANS.
  @param distance the distance functor.
  @param p the number of nearest neighbors of centers to estimate the width
  of Gaussian RBF functions.
  @return Gaussian RBF functions with parameter learned from data.