def METHOD_NAME(
    time_series: np.ndarray, minimum_samples: int = 3
) -> TimeSeriesStatistics:
    """Detect when a time series set has effectively become stationary (i.e has reached
    equilibrium).
    Parameters
    ----------
    time_series
        The time series to analyse with shape=(n_data_points, n_dimensions).
    minimum_samples
        The minimum number of data points to consider in the calculation.
    Notes
    -----
    This method is based on the paper by J. D. Chodera [1]_ and the implementation at
    https://github.com/choderalab/pymbar. Here the code is extended support
    multidimensional data such as dipole moments.
    References
    ----------
    [1] J. D. Chodera, W. C. Swope, J. W. Pitera, C. Seok, and K. A. Dill. Use of the
        weighted histogram analysis method for the analysis of simulated and parallel
        tempering simulations. JCTC 3(1):26-41, 2007.
    Returns
    -------
        Statistics about the time series.
    """
    n_timesteps = time_series.shape[0]
    statistical_inefficiency_array = np.ones([n_timesteps - 1])
    # Special case if the time series is constant.
    if np.isclose(time_series.std(), 0.0):
        return TimeSeriesStatistics(
            n_total_points=len(time_series),
            n_uncorrelated_points=1,
            statistical_inefficiency=float(len(time_series)),
            equilibration_index=0,
        )
    effect_samples_array = np.ones([n_timesteps - 1])
    for current_timestep in range(0, n_timesteps - 1):
        try:
            statistical_inefficiency_array[
                current_timestep
            ] = _statistical_inefficiency(
                time_series[current_timestep:n_timesteps], minimum_samples
            )
        except ParameterError:
            # Fix for issue https://github.com/choderalab/pymbar/issues/122
            statistical_inefficiency_array[current_timestep] = (
                n_timesteps - current_timestep + 1
            )
        effect_samples_array[current_timestep] = (
            n_timesteps - current_timestep + 1
        ) / statistical_inefficiency_array[current_timestep]
    equilibration_time = effect_samples_array.argmax()
    statistical_inefficiency = statistical_inefficiency_array[equilibration_time]
    return TimeSeriesStatistics(
        n_total_points=len(time_series),
        n_uncorrelated_points=len(
            get_uncorrelated_indices(len(time_series), statistical_inefficiency)
        ),
        statistical_inefficiency=float(statistical_inefficiency),
        equilibration_index=int(equilibration_time),
    )