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.

def choose_solution_using_percentiles(
        X_original,
        solutions,
        parameters=None,
        verbose=False,
        percentiles=list(range(10, 100, 10))):
    """
    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.
    """
    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