code: def choose_solution_using_percentiles(
          X_original,
          solutions,
          parameters=None,
          verbose=False,
          percentiles=list(range(10, 100, 10))):
      
      missing_mask = np.isnan(X_original)
      min_mse = np.inf
      best_solution = None
      for i, candidate in enumerate(solutions):
          for col_idx in range(X_original.shape[1]):
              col_data = candidate[:, col_idx]
              col_missing = missing_mask[:, col_idx]
              col_observed = ~col_missing
              if col_missing.sum() < 2:
                  continue
              elif col_observed.sum() < 2:
                  continue
              missing_data = col_data[col_missing]
              observed_data = col_data[col_observed]
  
              missing_percentiles = np.array([
                  np.percentile(missing_data, p)
                  for p in percentiles])
  
              observed_percentiles = np.array([
                  np.percentile(observed_data, p)
                  for p in percentiles])
  
              mse = np.mean((missing_percentiles - observed_percentiles) ** 2)
          if mse < min_mse:
              min_mse = mse
              best_solution = candidate
          if verbose:
              print("Candidate #%d/%d%s: %f" % (
                  i + 1,
                  len(solutions),
                  (" (parameter=%s) " % parameters[i]
                      if parameters is not None
                      else ""),
                  mse))
      return best_solution
docs: It's tricky to pick a single matrix out of all the candidate
      solutions with differing shrinkage thresholds.
      Our heuristic is to pick the matrix whose percentiles match best
      between the missing and observed data.