def estimate_multi_temperature(
    energy_trajs, temp_trajs, dtrajs,
    energy_unit='kcal/mol', temp_unit='K', reference_temperature=None,
    maxiter=10000, maxerr=1.0E-15, save_convergence_info=0,
    estimator='wham', lag=1, dt_traj='1 step', init=None, init_maxiter=10000, init_maxerr=1e-8,
    **kwargs):
    r"""
    This function acts as a wrapper for ``tram()``, ``dtram()``, ``mbar``, and ``wham()`` and
    handles the calculation of bias energies (``bias``) and thermodynamic state trajectories
    (``ttrajs``) when the data comes from multi-temperature simulations.

    Parameters
    ----------
    energy_trajs : list of N arrays, each of shape (T_i,)
        List of arrays, each having T_i rows, one for each time step, containing the potential
        energies time series in units of kT, kcal/mol or kJ/mol.
    temp_trajs : list of N int arrays, each of shape (T_i,)
        List of arrays, each having T_i rows, one for each time step, containing the heat bath
        temperature time series (at which temperatures the frames were created) in units of K or C.
        Alternatively, these trajectories may contain kT values instead of temperatures.
    dtrajs : list of N int arrays, each of shape (T_i,)
        The integers are indexes in 0,...,n-1 enumerating the n discrete states or the bins the
        trajectory is in at any time.
    energy_unit: str, optional, default='kcal/mol'
        The physical unit used for energies. Current options: kcal/mol, kJ/mol, kT.
    temp_unit : str, optional, default='K'
        The physical unit used for the temperature. Current options: K, C, kT
    reference_temperature : float or None, optional, default=None
        Reference temperature against which the bias energies are computed. If not given, the lowest
        temperature or kT value is used. If given, this parameter must have the same unit as the
        temp_trajs.
    maxiter : int, optional, default=10000
        The maximum number of self-consistent iterations before the estimator exits unsuccessfully.
    maxerr : float, optional, default=1E-15
        Convergence criterion based on the maximal free energy change in a self-consistent
        iteration step.
    save_convergence_info : int, optional, default=0
        Every save_convergence_info iteration steps, store the actual increment
        and the actual loglikelihood; 0 means no storage.
    estimator : str, optional, default='wham'
        Specify one of the available estimators

        | 'wham':   use WHAM
        | 'mbar':   use MBAR
        | 'dtram':  use the discrete version of TRAM
        | 'tram':  use TRAM
    lag : int or list of int, optional, default=1
        Integer lag time at which transitions are counted. Providing a list of lag times will
        trigger one estimation per lag time.
    dt_traj : str, optional, default='1 step'
        Description of the physical time corresponding to the lag. May be used by analysis
        algorithms such as plotting tools to pretty-print the axes. By default '1 step', i.e.
        there is no physical time unit.  Specify by a number, whitespace and unit. Permitted
        units are (* is an arbitrary string):

        |  'fs',   'femtosecond*'
        |  'ps',   'picosecond*'
        |  'ns',   'nanosecond*'
        |  'us',   'microsecond*'
        |  'ms',   'millisecond*'
        |  's',    'second*'
    init : str, optional, default=None
        Use a specific initialization for the self-consistent iteration:

        | None:    use a hard-coded guess for free energies and Lagrangian multipliers
        | 'wham':  perform a short WHAM estimate to initialize the free energies (only with dtram)
        | 'mbar':  perform a short MBAR estimate to initialize the free energies (only with tram)
    init_maxiter : int, optional, default=10000
        The maximum number of self-consistent iterations during the initialization.
    init_maxerr : float, optional, default=1.0E-8
        Convergence criterion for the initialization.
    **kwargs : dict, optional
        You can use this to pass estimator-specific named parameters to the chosen estimator, which
        are not already coverd by ``estimate_multi_temperature()``.

    Returns
    -------
    A :class:`MultiThermModel <pyemma.thermo.models.multi_therm.MultiThermModel>` or :class:`MEMM <pyemma.thermo.models.memm.MEMM>` object or list thereof
        The requested estimator/model object, i.e., WHAM, MBAR, DTRAM or TRAM. If multiple lag times
        are given, a list of objects is returned (one MEMM per lag time).

    Example
    -------
    We look at 1D simulations at two different kT values 1.0 and 2.0, already clustered data, and
    we use TRAM for the estimation:

    >>> from pyemma.thermo import estimate_multi_temperature as estimate_mt
    >>> import numpy as np
    >>> energy_trajs = [np.array([1.6, 1.4, 1.0, 1.0, 1.2, 1.0, 1.0]), np.array([0.8, 0.7, 0.5, 0.6, 0.7, 0.8, 0.7])]
    >>> temp_trajs = [np.array([1.0, 1.0, 1.0, 1.0, 1.0, 1.0, 1.0]), np.array([2.0, 2.0, 2.0, 2.0, 2.0, 2.0, 2.0])]
    >>> dtrajs = [np.array([0, 1, 2, 2, 2, 2, 2]), np.array([0, 1, 2, 2, 1, 0, 1])]
    >>> tram = estimate_mt(energy_trajs, temp_trajs, dtrajs, energy_unit='kT', temp_unit='kT', estimator='tram', lag=1)
    >>> tram.f # doctest: +ELLIPSIS
    array([ 2.90...,  1.72...,  0.26...])

    Note that alhough we only used one temperature per trajectory, ``estimate_multi_temperature()``
    can handle temperature changes as well.

    See :class:`MultiThermModel <pyemma.thermo.models.multi_therm.MultiThermModel>`
    or :class:`MEMM <pyemma.thermo.models.memm.MEMM>` for a full documentation.

    .. autoclass:: pyemma.thermo.models.multi_therm.MultiThermModel
        :members:
        :undoc-members:

        .. rubric:: Methods

        .. autoautosummary:: pyemma.thermo.models.multi_therm.MultiThermModel
           :methods:

        .. rubric:: Attributes

        .. autoautosummary:: pyemma.thermo.models.multi_therm.MultiThermModel
            :attributes:

    .. autoclass:: pyemma.thermo.models.memm.MEMM
        :members:
        :undoc-members:

        .. rubric:: Methods

        .. autoautosummary:: pyemma.thermo.models.memm.MEMM
           :methods:

        .. rubric:: Attributes

        .. autoautosummary:: pyemma.thermo.models.memm.MEMM
            :attributes:

    """
    if estimator not in ['wham', 'mbar', 'dtram', 'tram']:
        ValueError("unsupported estimator: %s" % estimator)
    from .util import get_multi_temperature_data as _get_multi_temperature_data
    ttrajs, btrajs, temperatures, unbiased_state = _get_multi_temperature_data(
        energy_trajs, temp_trajs, energy_unit, temp_unit,
        reference_temperature=reference_temperature)
    estimator_obj = None
    if estimator == 'wham':
        estimator_obj = wham(
            ttrajs, dtrajs,
            _get_averaged_bias_matrix(btrajs, dtrajs),
            maxiter=maxiter, maxerr=maxerr,
            save_convergence_info=save_convergence_info, dt_traj=dt_traj)
    elif estimator == 'mbar':
        allowed_keys = ['direct_space']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = mbar(
            ttrajs, dtrajs, btrajs,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, **parsed_kwargs)
    elif estimator == 'dtram':
        allowed_keys = ['count_mode', 'connectivity']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = dtram(
            ttrajs, dtrajs,
            _get_averaged_bias_matrix(btrajs, dtrajs),
            lag, unbiased_state=unbiased_state,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, init=init, init_maxiter=init_maxiter, init_maxerr=init_maxerr,
            **parsed_kwargs)
    elif estimator == 'tram':
        allowed_keys = [
            'count_mode', 'connectivity', 'connectivity_factor','nn',
            'direct_space', 'N_dtram_accelerations', 'equilibrium', 'overcounting_factor', 'callback']
        parsed_kwargs = dict([(i, kwargs[i]) for i in allowed_keys if i in kwargs])
        estimator_obj = tram(
            ttrajs, dtrajs, btrajs, lag, unbiased_state=unbiased_state,
            maxiter=maxiter, maxerr=maxerr, save_convergence_info=save_convergence_info,
            dt_traj=dt_traj, init=init, init_maxiter=init_maxiter, init_maxerr=init_maxerr,
            **parsed_kwargs)
    try:
        estimator_obj.temperatures = temperatures
    except AttributeError:
        for obj in estimator_obj:
            obj.temperatures = temperatures
    return estimator_obj