def select_levels(data, backtrack, kmin, kmax)
      return kmin if kmin == kmax

      method = :normal # "uniform" or "normal"

      kopt = kmin

      base = 1 # The position of first element in x: 1 or 0.
      n = data.size - base

      max_bic = 0.0

      for k in kmin..kmax
        cluster_sizes = []
        kbacktrack = backtrack[0..k]
        backtrack(kbacktrack) do |cluster, left, right|
          cluster_sizes[cluster] = right - left + 1
        end

        index_left = base
        index_right = 0

        likelihood = 0
        bin_left, bin_right = 0
        for i in 0..(k-1)
          points_in_bin = cluster_sizes[i + base]
          index_right = index_left + points_in_bin - 1

          if data[index_left] < data[index_right]
            bin_left = data[index_left]
            bin_right = data[index_right]
          elsif data[index_left] == data[index_right]
            bin_left = index_left == base ? data[base] : (data[index_left-1] + data[index_left]) / 2
            bin_right = index_right < n-1+base ? (data[index_right] + data[index_right+1]) / 2 : data[n-1+base]
          else
            raise "ERROR: binLeft > binRight"
          end

          bin_width = bin_right - bin_left
          if method == :uniform
            likelihood += points_in_bin * Math.log(points_in_bin / bin_width / n)
          else
            mean = 0.0
            variance = 0.0

            for j in index_left..index_right
              mean += data[j]
              variance += data[j] ** 2
            end
            mean /= points_in_bin
            variance = (variance - points_in_bin * mean ** 2) / (points_in_bin - 1) if points_in_bin > 1

            if variance > 0
              for j in index_left..index_right
                likelihood += - (data[j] - mean) ** 2 / (2.0 * variance)
              end
              likelihood += points_in_bin * (Math.log(points_in_bin / Float(n))
                - 0.5 * Math.log( 2 * Math::PI * variance))
            else
              likelihood += points_in_bin * Math.log(1.0 / bin_width / n)
            end
          end

          index_left = index_right + 1
        end

        # Compute the Bayesian information criterion
        bic = 2 * likelihood - (3 * k - 1) * Math.log(Float(n))

        if k == kmin
          max_bic = bic
          kopt = kmin
        elsif bic > max_bic
          max_bic = bic
          kopt = k
        end

      end

      kopt
    end
def select_levels(data, backtrack, kmin, kmax)
      return kmin if kmin == kmax

      method = :normal # "uniform" or "normal"

      kopt = kmin

      base = 1 # The position of first element in x: 1 or 0.
      n = data.size - base

      max_bic = 0.0

      for k in kmin..kmax
        cluster_sizes = []
        kbacktrack = backtrack[0..k]
        backtrack(kbacktrack) do |cluster, left, right|
          cluster_sizes[cluster] = right - left + 1
        end

        index_left = base
        index_right = 0

        likelihood = 0
        bin_left, bin_right = 0
        for i in 0..(k-1)
          points_in_bin = cluster_sizes[i + base]
          index_right = index_left + points_in_bin - 1

          if data[index_left] < data[index_right]
            bin_left = data[index_left]
            bin_right = data[index_right]
          elsif data[index_left] == data[index_right]
            bin_left = index_left == base ? data[base] : (data[index_left-1] + data[index_left]) / 2
            bin_right = index_right < n-1+base ? (data[index_right] + data[index_right+1]) / 2 : data[n-1+base]
          else
            raise "ERROR: binLeft > binRight"
          end

          bin_width = bin_right - bin_left
          if method == :uniform
            likelihood += points_in_bin * Math.log(points_in_bin / bin_width / n)
          else
            mean = 0.0
            variance = 0.0

            for j in index_left..index_right
              mean += data[j]
              variance += data[j] ** 2
            end
            mean /= points_in_bin
            variance = (variance - points_in_bin * mean ** 2) / (points_in_bin - 1) if points_in_bin > 1

            if variance > 0
              for j in index_left..index_right
                likelihood += - (data[j] - mean) ** 2 / (2.0 * variance)
              end
              likelihood += points_in_bin * (Math.log(points_in_bin / Float(n))
                - 0.5 * Math.log( 2 * Math::PI * variance))
            else
              likelihood += points_in_bin * Math.log(1.0 / bin_width / n)
            end
          end

          index_left = index_right + 1
        end

        # Compute the Bayesian information criterion
        bic = 2 * likelihood - (3 * k - 1) * Math.log(Float(n))

        if k == kmin
          max_bic = bic
          kopt = kmin
        elsif bic > max_bic
          max_bic = bic
          kopt = k
        end

      end

      kopt
    end
Choose an optimal number of levels between Kmin and Kmax