Repository: scikit-learn/scikit-learn
Issue category: Feature Request

Issue:
MICE imputer
Proceeding from https://github.com/scikit-learn/scikit-learn/pull/4844#issuecomment-253043238 is a suggestion that MICE imputation be included in scikit-learn. @sergeyf has implemented it [here](https://github.com/hammerlab/fancyimpute/blob/master/fancyimpute/mice.py).

Here we will discuss issues related to its design and inclusion.




Hints:
@sergeyf [wrote](https://github.com/scikit-learn/scikit-learn/pull/4844#issuecomment-259055975):

> @amueller and @jnothman do either of you have a sense of what the right structure for MICE would be? It's pretty involved, so I don't want to just add a new option to `Imputer`. What about adding a new `MICE` class to `imputation.py`?
> 
> Also, the `fit` method doesn't really make sense for MICE. It can pretty much only do `fit_transform` or `transform`: you provide some `X_train` that has holes in it, and MICE fills them in. If you then have an `X_test`, it doesn't use any of the info from `X_train` to fill it in. The simplest thing to do here would be to have `fit` do nothing, have `transform` do all the work, and have `fit_transform` alias to `transform`.

Yes, please create a new `MICEImputer` class.

I think we need an API that allows people to use this in a standard supervised prediction pipeline. Two issues with your proposal then:
1. the data being transformed is often much smaller (and therefore less informative) than the training data. To solve this, we could store the training data (pre-imputed, IMO) and perform the imputation rounds on the concatenation of training and transform data. However:
2. a downstream estimator will operate best when the training transformation and testing transformation is identical. To solve this, we could instead store the model constructed for each column in each of the `n_imputations` performed at train time. (In the case of PMM, we would draw the K nearest neighbors only from the training data as imputed at that point in time.) This also ensures that each row is transformed deterministically given a fitted model.

Could you please comment on how this relates to what's possible in the R interface?

@jnothman It is certainly often true that the data being transformed is smaller than the training data. I like the idea of storing each model for each iteration that was done during the fit phase (as well as the order the imputations were done in). So this would be `fit` or `fit_transform`?

Also, I wasn't going to implement PMM for this pull request: it's much slower and did not work nearly as well in most of my experiments. If someone requests it later, I can add it on without too much trouble.  

As for MICE in R: I'm not too familiar with it, but I think MICE just fills in matrices. There is no notion of training and test - that is left up to the user. 

One key point is that MICE provides the last `n_last` imputations, and doesn't average them for you. The intent is to then take the multiple imputations and use them all in parallel downstream, averaging only as the final step (after, say, training and prediction of `n_last` classifiers). We could return both an averaged version, and optionally return the last `n_last` imputations. What do you think?

`fit` and `fit_transform` are meant to be functionally identical in terms of the model that is stored. Both would run the full imputation process on the training data and store the sufficient information to run the imputations on new data in `transform`. I suspect in this case little would differ between their implementations but the return value.

OK, sounds good, thanks. 

I think I have a sense of what to do. I'll get started once the `return_std` is accepted.

Though you're right, we don't _need_ any ability to transform new data (just as we do not provide this in TSNE), as the technique is unsupervised. So first implement `fit_transform`, then we'll decide whether it's worth storing the model and implementing `transform`.

Well I already implemented `transform` as it was fairly straight-forward. I'll leave it in and we can trim it later if that seems like the wiser path. 


Patch:
diff --git a/doc/modules/classes.rst b/doc/modules/classes.rst
index 2e7dcba82e846..243c63ab0c7e2 100644
--- a/doc/modules/classes.rst
+++ b/doc/modules/classes.rst
@@ -646,6 +646,7 @@ Kernels:
    :template: class.rst
 
    impute.SimpleImputer
+   impute.MICEImputer
 
 .. _kernel_approximation_ref:
 
diff --git a/doc/modules/impute.rst b/doc/modules/impute.rst
index e806cc2fd5b4a..c4e8a3395025c 100644
--- a/doc/modules/impute.rst
+++ b/doc/modules/impute.rst
@@ -15,6 +15,10 @@ values. However, this comes at the price of losing data which may be valuable
 (even though incomplete). A better strategy is to impute the missing values,
 i.e., to infer them from the known part of the data.
 
+
+Univariate feature imputation
+=============================
+
 The :class:`SimpleImputer` class provides basic strategies for imputing missing
 values, either using the mean, the median or the most frequent value of
 the row or column in which the missing values are located. This class
@@ -52,5 +56,34 @@ Note that, here, missing values are encoded by 0 and are thus implicitly stored
 in the matrix. This format is thus suitable when there are many more missing
 values than observed values.
 
-:class:`SimpleImputer` can be used in a Pipeline as a way to build a composite
-estimator that supports imputation. See :ref:`sphx_glr_auto_examples_plot_missing_values.py`.
+
+Multivariate feature imputation
+===============================
+
+A more sophisticated approach is to use the :class:`MICEImputer` class, which
+implements the Multivariate Imputation by Chained Equations technique. MICE
+models each feature with missing values as a function of other features, and
+uses that estimate for imputation. It does so in a round-robin fashion: at
+each step, a feature column is designated as output `y` and the other feature
+columns are treated as inputs `X`. A regressor is fit on `(X, y)` for known `y`.
+Then, the regressor is used to predict the unknown values of `y`. This is
+repeated for each feature, and then is done for a number of imputation rounds.
+Here is an example snippet::
+
+    >>> import numpy as np
+    >>> from sklearn.impute import MICEImputer
+    >>> imp = MICEImputer(n_imputations=10, random_state=0)
+    >>> imp.fit([[1, 2], [np.nan, 3], [7, np.nan]])
+    MICEImputer(imputation_order='ascending', initial_strategy='mean',
+          max_value=None, min_value=None, missing_values='NaN', n_burn_in=10,
+          n_imputations=10, n_nearest_features=None, predictor=None,
+          random_state=0, verbose=False)
+    >>> X_test = [[np.nan, 2], [6, np.nan], [np.nan, 6]]
+    >>> print(np.round(imp.transform(X_test)))
+    [[ 1.  2.]
+     [ 6.  4.]
+     [13.  6.]]
+
+Both :class:`SimpleImputer` and :class:`MICEImputer` can be used in a Pipeline
+as a way to build a composite estimator that supports imputation.
+See :ref:`sphx_glr_auto_examples_plot_missing_values.py`.
diff --git a/examples/plot_missing_values.py b/examples/plot_missing_values.py
index 9cc07b75d2542..e32c19fae0846 100644
--- a/examples/plot_missing_values.py
+++ b/examples/plot_missing_values.py
@@ -1,72 +1,127 @@
 """
-======================================================
+====================================================
 Imputing missing values before building an estimator
-======================================================
-
-This example shows that imputing the missing values can give better
-results than discarding the samples containing any missing value.
-Imputing does not always improve the predictions, so please check via
-cross-validation.  Sometimes dropping rows or using marker values is
-more effective.
+====================================================
 
 Missing values can be replaced by the mean, the median or the most frequent
-value using the ``strategy`` hyper-parameter.
+value using the basic ``SimpleImputer``.
 The median is a more robust estimator for data with high magnitude variables
 which could dominate results (otherwise known as a 'long tail').
 
-Script output::
-
-  Score with the entire dataset = 0.56
-  Score without the samples containing missing values = 0.48
-  Score after imputation of the missing values = 0.55
-
-In this case, imputing helps the classifier get close to the original score.
-
+Another option is the MICE imputer. This uses round-robin linear regression,
+treating every variable as an output in turn. The version implemented assumes
+Gaussian (output) variables. If your features are obviously non-Normal,
+consider transforming them to look more Normal so as to improve performance.
 """
+
 import numpy as np
+import matplotlib.pyplot as plt
 
+from sklearn.datasets import load_diabetes
 from sklearn.datasets import load_boston
 from sklearn.ensemble import RandomForestRegressor
 from sklearn.pipeline import Pipeline
-from sklearn.impute import SimpleImputer
+from sklearn.impute import SimpleImputer, MICEImputer
 from sklearn.model_selection import cross_val_score
 
 rng = np.random.RandomState(0)
 
-dataset = load_boston()
-X_full, y_full = dataset.data, dataset.target
-n_samples = X_full.shape[0]
-n_features = X_full.shape[1]
-
-# Estimate the score on the entire dataset, with no missing values
-estimator = RandomForestRegressor(random_state=0, n_estimators=100)
-score = cross_val_score(estimator, X_full, y_full).mean()
-print("Score with the entire dataset = %.2f" % score)
-
-# Add missing values in 75% of the lines
-missing_rate = 0.75
-n_missing_samples = int(np.floor(n_samples * missing_rate))
-missing_samples = np.hstack((np.zeros(n_samples - n_missing_samples,
-                                      dtype=np.bool),
-                             np.ones(n_missing_samples,
-                                     dtype=np.bool)))
-rng.shuffle(missing_samples)
-missing_features = rng.randint(0, n_features, n_missing_samples)
-
-# Estimate the score without the lines containing missing values
-X_filtered = X_full[~missing_samples, :]
-y_filtered = y_full[~missing_samples]
-estimator = RandomForestRegressor(random_state=0, n_estimators=100)
-score = cross_val_score(estimator, X_filtered, y_filtered).mean()
-print("Score without the samples containing missing values = %.2f" % score)
-
-# Estimate the score after imputation of the missing values
-X_missing = X_full.copy()
-X_missing[np.where(missing_samples)[0], missing_features] = 0
-y_missing = y_full.copy()
-estimator = Pipeline([("imputer", SimpleImputer(missing_values=0,
-                                                strategy="mean")),
-                      ("forest", RandomForestRegressor(random_state=0,
-                                                       n_estimators=100))])
-score = cross_val_score(estimator, X_missing, y_missing).mean()
-print("Score after imputation of the missing values = %.2f" % score)
+
+def get_results(dataset):
+    X_full, y_full = dataset.data, dataset.target
+    n_samples = X_full.shape[0]
+    n_features = X_full.shape[1]
+
+    # Estimate the score on the entire dataset, with no missing values
+    estimator = RandomForestRegressor(random_state=0, n_estimators=100)
+    full_scores = cross_val_score(estimator, X_full, y_full,
+                                  scoring='neg_mean_squared_error')
+
+    # Add missing values in 75% of the lines
+    missing_rate = 0.75
+    n_missing_samples = int(np.floor(n_samples * missing_rate))
+    missing_samples = np.hstack((np.zeros(n_samples - n_missing_samples,
+                                          dtype=np.bool),
+                                 np.ones(n_missing_samples,
+                                         dtype=np.bool)))
+    rng.shuffle(missing_samples)
+    missing_features = rng.randint(0, n_features, n_missing_samples)
+
+    # Estimate the score after replacing missing values by 0
+    X_missing = X_full.copy()
+    X_missing[np.where(missing_samples)[0], missing_features] = 0
+    y_missing = y_full.copy()
+    estimator = RandomForestRegressor(random_state=0, n_estimators=100)
+    zero_impute_scores = cross_val_score(estimator, X_missing, y_missing,
+                                         scoring='neg_mean_squared_error')
+
+    # Estimate the score after imputation (mean strategy) of the missing values
+    X_missing = X_full.copy()
+    X_missing[np.where(missing_samples)[0], missing_features] = 0
+    y_missing = y_full.copy()
+    estimator = Pipeline([("imputer", SimpleImputer(missing_values=0,
+                                                    strategy="mean")),
+                          ("forest", RandomForestRegressor(random_state=0,
+                                                           n_estimators=100))])
+    mean_impute_scores = cross_val_score(estimator, X_missing, y_missing,
+                                         scoring='neg_mean_squared_error')
+
+    # Estimate the score after imputation (MICE strategy) of the missing values
+    estimator = Pipeline([("imputer", MICEImputer(missing_values=0,
+                                                  random_state=0)),
+                          ("forest", RandomForestRegressor(random_state=0,
+                                                           n_estimators=100))])
+    mice_impute_scores = cross_val_score(estimator, X_missing, y_missing,
+                                         scoring='neg_mean_squared_error')
+
+    return ((full_scores.mean(), full_scores.std()),
+            (zero_impute_scores.mean(), zero_impute_scores.std()),
+            (mean_impute_scores.mean(), mean_impute_scores.std()),
+            (mice_impute_scores.mean(), mice_impute_scores.std()))
+
+
+results_diabetes = np.array(get_results(load_diabetes()))
+mses_diabetes = results_diabetes[:, 0] * -1
+stds_diabetes = results_diabetes[:, 1]
+
+results_boston = np.array(get_results(load_boston()))
+mses_boston = results_boston[:, 0] * -1
+stds_boston = results_boston[:, 1]
+
+n_bars = len(mses_diabetes)
+xval = np.arange(n_bars)
+
+x_labels = ['Full data',
+            'Zero imputation',
+            'Mean Imputation',
+            'MICE Imputation']
+colors = ['r', 'g', 'b', 'orange']
+
+# plot diabetes results
+plt.figure(figsize=(12, 6))
+ax1 = plt.subplot(121)
+for j in xval:
+    ax1.barh(j, mses_diabetes[j], xerr=stds_diabetes[j],
+             color=colors[j], alpha=0.6, align='center')
+
+ax1.set_title('Feature Selection Techniques with Diabetes Data')
+ax1.set_xlim(left=np.min(mses_diabetes) * 0.9,
+             right=np.max(mses_diabetes) * 1.1)
+ax1.set_yticks(xval)
+ax1.set_xlabel('MSE')
+ax1.invert_yaxis()
+ax1.set_yticklabels(x_labels)
+
+# plot boston results
+ax2 = plt.subplot(122)
+for j in xval:
+    ax2.barh(j, mses_boston[j], xerr=stds_boston[j],
+             color=colors[j], alpha=0.6, align='center')
+
+ax2.set_title('Feature Selection Techniques with Boston Data')
+ax2.set_yticks(xval)
+ax2.set_xlabel('MSE')
+ax2.invert_yaxis()
+ax2.set_yticklabels([''] * n_bars)
+
+plt.show()
diff --git a/sklearn/impute.py b/sklearn/impute.py
index 2420e02560e42..f395c634ccced 100644
--- a/sklearn/impute.py
+++ b/sklearn/impute.py
@@ -1,16 +1,23 @@
 """Transformers for missing value imputation"""
 # Authors: Nicolas Tresegnie <nicolas.tresegnie@gmail.com>
+#          Sergey Feldman <sergeyfeldman@gmail.com>
 # License: BSD 3 clause
 
+from __future__ import division
+
 import warnings
+from time import time
 
 import numpy as np
 import numpy.ma as ma
 from scipy import sparse
 from scipy import stats
+from collections import namedtuple
 
 from .base import BaseEstimator, TransformerMixin
-from .utils import check_array
+from .base import clone
+from .preprocessing import normalize
+from .utils import check_array, check_random_state, safe_indexing
 from .utils.sparsefuncs import _get_median
 from .utils.validation import check_is_fitted
 from .utils.validation import FLOAT_DTYPES
@@ -20,8 +27,13 @@
 zip = six.moves.zip
 map = six.moves.map
 
+MICETriplet = namedtuple('MICETriplet', ['feat_idx',
+                                         'neighbor_feat_idx',
+                                         'predictor'])
+
 __all__ = [
     'SimpleImputer',
+    'MICEImputer',
 ]
 
 
@@ -321,3 +333,541 @@ def transform(self, X):
             X[coordinates] = values
 
         return X
+
+
+class MICEImputer(BaseEstimator, TransformerMixin):
+    """MICE transformer to impute missing values.
+
+    Basic implementation of MICE (Multivariate Imputations by Chained
+    Equations) package from R. This version assumes all of the features are
+    Gaussian.
+
+    Parameters
+    ----------
+    missing_values : int or "NaN", optional (default="NaN")
+        The placeholder for the missing values. All occurrences of
+        ``missing_values`` will be imputed. For missing values encoded as
+        np.nan, use the string value "NaN".
+
+    imputation_order : str, optional (default="ascending")
+        The order in which the features will be imputed. Possible values:
+
+        "ascending"
+            From features with fewest missing values to most.
+        "descending"
+            From features with most missing values to fewest.
+        "roman"
+            Left to right.
+        "arabic"
+            Right to left.
+        "random"
+            A random order for each round.
+
+    n_imputations : int, optional (default=100)
+        Number of MICE rounds to perform, the results of which will be
+        used in the final average.
+
+    n_burn_in : int, optional (default=10)
+        Number of initial MICE rounds to perform the results of which
+        will not be returned.
+
+    predictor : estimator object, default=BayesianRidge()
+        The predictor to use at each step of the round-robin imputation.
+        It must support ``return_std`` in its ``predict`` method.
+
+    n_nearest_features : int, optional (default=None)
+        Number of other features to use to estimate the missing values of
+        the each feature column. Nearness between features is measured using
+        the absolute correlation coefficient between each feature pair (after
+        initial imputation). Can provide significant speed-up when the number
+        of features is huge. If ``None``, all features will be used.
+
+    initial_strategy : str, optional (default="mean")
+        Which strategy to use to initialize the missing values. Same as the
+        ``strategy`` parameter in :class:`sklearn.preprocessing.Imputer`
+        Valid values: {"mean", "median", or "most_frequent"}.
+
+    min_value : float, optional (default=None)
+        Minimum possible imputed value. Default of ``None`` will set minimum
+        to negative infinity.
+
+    max_value : float, optional (default=None)
+        Maximum possible imputed value. Default of ``None`` will set maximum
+        to positive infinity.
+
+    verbose : int, optional (default=0)
+        Verbosity flag, controls the debug messages that are issued
+        as functions are evaluated. The higher, the more verbose. Can be 0, 1,
+        or 2.
+
+    random_state : int, RandomState instance or None, optional (default=None)
+        The seed of the pseudo random number generator to use when shuffling
+        the data.  If int, random_state is the seed used by the random number
+        generator; If RandomState instance, random_state is the random number
+        generator; If None, the random number generator is the RandomState
+        instance used by ``np.random``.
+
+    Attributes
+    ----------
+    initial_imputer_ : object of class :class:`sklearn.preprocessing.Imputer`'
+        The imputer used to initialize the missing values.
+
+    imputation_sequence_ : list of tuples
+        Each tuple has ``(feat_idx, neighbor_feat_idx, predictor)``, where
+        ``feat_idx`` is the current feature to be imputed,
+        ``neighbor_feat_idx`` is the array of other features used to impute the
+        current feature, and ``predictor`` is the trained predictor used for
+        the imputation.
+
+    Notes
+    -----
+    The R version of MICE does not have inductive functionality, i.e. first
+    fitting on ``X_train`` and then transforming any ``X_test`` without
+    additional fitting. We do this by storing each feature's predictor during
+    the round-robin ``fit`` phase, and predicting without refitting (in order)
+    during the ``transform`` phase.
+
+    Features which contain all missing values at ``fit`` are discarded upon
+    ``transform``.
+
+    Features with missing values in transform which did not have any missing
+    values in fit will be imputed with the initial imputation method only.
+
+    References
+    ----------
+    .. [1] `Stef van Buuren, Karin Groothuis-Oudshoorn (2011). "mice:
+        Multivariate Imputation by Chained Equations in R". Journal of
+        Statistical Software 45: 1-67.
+        <https://www.jstatsoft.org/article/view/v045i03>`_
+    """
+
+    def __init__(self,
+                 missing_values='NaN',
+                 imputation_order='ascending',
+                 n_imputations=100,
+                 n_burn_in=10,
+                 predictor=None,
+                 n_nearest_features=None,
+                 initial_strategy="mean",
+                 min_value=None,
+                 max_value=None,
+                 verbose=False,
+                 random_state=None):
+
+        self.missing_values = missing_values
+        self.imputation_order = imputation_order
+        self.n_imputations = n_imputations
+        self.n_burn_in = n_burn_in
+        self.predictor = predictor
+        self.n_nearest_features = n_nearest_features
+        self.initial_strategy = initial_strategy
+        self.min_value = min_value
+        self.max_value = max_value
+        self.verbose = verbose
+        self.random_state = random_state
+
+    def _impute_one_feature(self,
+                            X_filled,
+                            mask_missing_values,
+                            feat_idx,
+                            neighbor_feat_idx,
+                            predictor=None,
+                            fit_mode=True):
+        """Impute a single feature from the others provided.
+
+        This function predicts the missing values of one of the features using
+        the current estimates of all the other features. The ``predictor`` must
+        support ``return_std=True`` in its ``predict`` method for this function
+        to work.
+
+        Parameters
+        ----------
+        X_filled : ndarray
+            Input data with the most recent imputations.
+
+        mask_missing_values : ndarray
+            Input data's missing indicator matrix.
+
+        feat_idx : int
+            Index of the feature currently being imputed.
+
+        neighbor_feat_idx : ndarray
+            Indices of the features to be used in imputing ``feat_idx``.
+
+        predictor : object
+            The predictor to use at this step of the round-robin imputation.
+            It must support ``return_std`` in its ``predict`` method.
+            If None, it will be cloned from self._predictor.
+
+        fit_mode : boolean, default=True
+            Whether to fit and predict with the predictor or just predict.
+
+        Returns
+        -------
+        X_filled : ndarray
+            Input data with ``X_filled[missing_row_mask, feat_idx]`` updated.
+
+        predictor : predictor with sklearn API
+            The fitted predictor used to impute
+            ``X_filled[missing_row_mask, feat_idx]``.
+        """
+
+        # if nothing is missing, just return the default
+        # (should not happen at fit time because feat_ids would be excluded)
+        missing_row_mask = mask_missing_values[:, feat_idx]
+        if not np.any(missing_row_mask):
+            return X_filled, predictor
+
+        if predictor is None and fit_mode is False:
+            raise ValueError("If fit_mode is False, then an already-fitted "
+                             "predictor should be passed in.")
+
+        if predictor is None:
+            predictor = clone(self._predictor)
+
+        if fit_mode:
+            X_train = safe_indexing(X_filled[:, neighbor_feat_idx],
+                                    ~missing_row_mask)
+            y_train = safe_indexing(X_filled[:, feat_idx],
+                                    ~missing_row_mask)
+            predictor.fit(X_train, y_train)
+
+        # get posterior samples
+        X_test = safe_indexing(X_filled[:, neighbor_feat_idx],
+                               missing_row_mask)
+        mus, sigmas = predictor.predict(X_test, return_std=True)
+        good_sigmas = sigmas > 0
+        imputed_values = np.zeros(mus.shape, dtype=X_filled.dtype)
+        imputed_values[~good_sigmas] = mus[~good_sigmas]
+        imputed_values[good_sigmas] = self.random_state_.normal(
+            loc=mus[good_sigmas], scale=sigmas[good_sigmas])
+
+        # clip the values
+        imputed_values = np.clip(imputed_values,
+                                 self._min_value,
+                                 self._max_value)
+
+        # update the feature
+        X_filled[missing_row_mask, feat_idx] = imputed_values
+        return X_filled, predictor
+
+    def _get_neighbor_feat_idx(self,
+                               n_features,
+                               feat_idx,
+                               abs_corr_mat):
+        """Get a list of other features to predict ``feat_idx``.
+
+        If self.n_nearest_features is less than or equal to the total
+        number of features, then use a probability proportional to the absolute
+        correlation between ``feat_idx`` and each other feature to randomly
+        choose a subsample of the other features (without replacement).
+
+        Parameters
+        ----------
+        n_features : int
+            Number of features in ``X``.
+
+        feat_idx : int
+            Index of the feature currently being imputed.
+
+        abs_corr_mat : ndarray, shape (n_features, n_features)
+            Absolute correlation matrix of ``X``. The diagonal has been zeroed
+            out and each feature has been normalized to sum to 1. Can be None.
+
+        Returns
+        -------
+        neighbor_feat_idx : array-like
+            The features to use to impute ``feat_idx``.
+        """
+        if (self.n_nearest_features is not None and
+                self.n_nearest_features < n_features):
+            p = abs_corr_mat[:, feat_idx]
+            neighbor_feat_idx = self.random_state_.choice(
+                np.arange(n_features), self.n_nearest_features, replace=False,
+                p=p)
+        else:
+            inds_left = np.arange(feat_idx)
+            inds_right = np.arange(feat_idx + 1, n_features)
+            neighbor_feat_idx = np.concatenate((inds_left, inds_right))
+        return neighbor_feat_idx
+
+    def _get_ordered_idx(self, mask_missing_values):
+        """Decide in what order we will update the features.
+
+        As a homage to the MICE R package, we will have 4 main options of
+        how to order the updates, and use a random order if anything else
+        is specified.
+
+        Also, this function skips features which have no missing values.
+
+        Parameters
+        ----------
+        mask_missing_values : array-like, shape (n_samples, n_features)
+            Input data's missing indicator matrix, where "n_samples" is the
+            number of samples and "n_features" is the number of features.
+
+        Returns
+        -------
+        ordered_idx : ndarray, shape (n_features,)
+            The order in which to impute the features.
+        """
+        frac_of_missing_values = mask_missing_values.mean(axis=0)
+        missing_values_idx = np.nonzero(frac_of_missing_values)[0]
+        if self.imputation_order == 'roman':
+            ordered_idx = missing_values_idx
+        elif self.imputation_order == 'arabic':
+            ordered_idx = missing_values_idx[::-1]
+        elif self.imputation_order == 'ascending':
+            n = len(frac_of_missing_values) - len(missing_values_idx)
+            ordered_idx = np.argsort(frac_of_missing_values,
+                                     kind='mergesort')[n:][::-1]
+        elif self.imputation_order == 'descending':
+            n = len(frac_of_missing_values) - len(missing_values_idx)
+            ordered_idx = np.argsort(frac_of_missing_values,
+                                     kind='mergesort')[n:]
+        elif self.imputation_order == 'random':
+            ordered_idx = missing_values_idx
+            self.random_state_.shuffle(ordered_idx)
+        else:
+            raise ValueError("Got an invalid imputation order: '{0}'. It must "
+                             "be one of the following: 'roman', 'arabic', "
+                             "'ascending', 'descending', or "
+                             "'random'.".format(self.imputation_order))
+        return ordered_idx
+
+    def _get_abs_corr_mat(self, X_filled, tolerance=1e-6):
+        """Get absolute correlation matrix between features.
+
+        Parameters
+        ----------
+        X_filled : ndarray, shape (n_samples, n_features)
+            Input data with the most recent imputations.
+
+        tolerance : float, optional (default=1e-6)
+            ``abs_corr_mat`` can have nans, which will be replaced
+            with ``tolerance``.
+
+        Returns
+        -------
+        abs_corr_mat : ndarray, shape (n_features, n_features)
+            Absolute correlation matrix of ``X`` at the beginning of the
+            current round. The diagonal has been zeroed out and each feature's
+            absolute correlations with all others have been normalized to sum
+            to 1.
+        """
+        n_features = X_filled.shape[1]
+        if (self.n_nearest_features is None or
+                self.n_nearest_features >= n_features):
+            return None
+        abs_corr_mat = np.abs(np.corrcoef(X_filled.T))
+        # np.corrcoef is not defined for features with zero std
+        abs_corr_mat[np.isnan(abs_corr_mat)] = tolerance
+        # ensures exploration, i.e. at least some probability of sampling
+        abs_corr_mat[abs_corr_mat < tolerance] = tolerance
+        # features are not their own neighbors
+        np.fill_diagonal(abs_corr_mat, 0)
+        # needs to sum to 1 for np.random.choice sampling
+        abs_corr_mat = normalize(abs_corr_mat, norm='l1', axis=0, copy=False)
+        return abs_corr_mat
+
+    def _initial_imputation(self, X):
+        """Perform initial imputation for input X.
+
+        Parameters
+        ----------
+        X : ndarray, shape (n_samples, n_features)
+            Input data, where "n_samples" is the number of samples and
+            "n_features" is the number of features.
+
+        Returns
+        -------
+        Xt : ndarray, shape (n_samples, n_features)
+            Input data, where "n_samples" is the number of samples and
+            "n_features" is the number of features.
+
+        X_filled : ndarray, shape (n_samples, n_features)
+            Input data with the most recent imputations.
+
+        mask_missing_values : ndarray, shape (n_samples, n_features)
+            Input data's missing indicator matrix, where "n_samples" is the
+            number of samples and "n_features" is the number of features.
+        """
+        X = check_array(X, dtype=FLOAT_DTYPES, order="F",
+                        force_all_finite='allow-nan'
+                        if self.missing_values == 'NaN'
+                        or np.isnan(self.missing_values) else True)
+
+        mask_missing_values = _get_mask(X, self.missing_values)
+        if self.initial_imputer_ is None:
+            self.initial_imputer_ = SimpleImputer(
+                                            missing_values=self.missing_values,
+                                            strategy=self.initial_strategy)
+            X_filled = self.initial_imputer_.fit_transform(X)
+        else:
+            X_filled = self.initial_imputer_.transform(X)
+
+        valid_mask = np.flatnonzero(np.logical_not(
+            np.isnan(self.initial_imputer_.statistics_)))
+        Xt = X[:, valid_mask]
+        mask_missing_values = mask_missing_values[:, valid_mask]
+
+        return Xt, X_filled, mask_missing_values
+
+    def fit_transform(self, X, y=None):
+        """Fits the imputer on X and return the transformed X.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            Input data, where "n_samples" is the number of samples and
+            "n_features" is the number of features.
+
+        y : ignored.
+
+        Returns
+        -------
+        Xt : array-like, shape (n_samples, n_features)
+             The imputed input data.
+        """
+        self.random_state_ = getattr(self, "random_state_",
+                                     check_random_state(self.random_state))
+
+        if self.predictor is None:
+            from .linear_model import BayesianRidge
+            self._predictor = BayesianRidge()
+        else:
+            self._predictor = clone(self.predictor)
+
+        self._min_value = np.nan if self.min_value is None else self.min_value
+        self._max_value = np.nan if self.max_value is None else self.max_value
+
+        self.initial_imputer_ = None
+        X, X_filled, mask_missing_values = self._initial_imputation(X)
+
+        # edge case: in case the user specifies 0 for n_imputations,
+        # then there is no need to do burn in and the result should be
+        # just the initial imputation (before clipping)
+        if self.n_imputations < 1:
+            return X_filled
+
+        X_filled = np.clip(X_filled, self._min_value, self._max_value)
+
+        # order in which to impute
+        # note this is probably too slow for large feature data (d > 100000)
+        # and a better way would be good.
+        # see: https://goo.gl/KyCNwj and subsequent comments
+        ordered_idx = self._get_ordered_idx(mask_missing_values)
+
+        abs_corr_mat = self._get_abs_corr_mat(X_filled)
+
+        # impute data
+        n_rounds = self.n_burn_in + self.n_imputations
+        n_samples, n_features = X_filled.shape
+        Xt = np.zeros((n_samples, n_features), dtype=X.dtype)
+        self.imputation_sequence_ = []
+        if self.verbose > 0:
+            print("[MICE] Completing matrix with shape %s" % (X.shape,))
+        start_t = time()
+        for i_rnd in range(n_rounds):
+            if self.imputation_order == 'random':
+                ordered_idx = self._get_ordered_idx(mask_missing_values)
+
+            for feat_idx in ordered_idx:
+                neighbor_feat_idx = self._get_neighbor_feat_idx(n_features,
+                                                                feat_idx,
+                                                                abs_corr_mat)
+                X_filled, predictor = self._impute_one_feature(
+                    X_filled, mask_missing_values, feat_idx, neighbor_feat_idx,
+                    predictor=None, fit_mode=True)
+                predictor_triplet = MICETriplet(feat_idx,
+                                                neighbor_feat_idx,
+                                                predictor)
+                self.imputation_sequence_.append(predictor_triplet)
+
+            if i_rnd >= self.n_burn_in:
+                Xt += X_filled
+            if self.verbose > 0:
+                print('[MICE] Ending imputation round '
+                      '%d/%d, elapsed time %0.2f'
+                      % (i_rnd + 1, n_rounds, time() - start_t))
+
+        Xt /= self.n_imputations
+        Xt[~mask_missing_values] = X[~mask_missing_values]
+        return Xt
+
+    def transform(self, X):
+        """Imputes all missing values in X.
+
+        Note that this is stochastic, and that if random_state is not fixed,
+        repeated calls, or permuted input, will yield different results.
+
+        Parameters
+        ----------
+        X : array-like, shape = [n_samples, n_features]
+            The input data to complete.
+
+        Returns
+        -------
+        Xt : array-like, shape (n_samples, n_features)
+             The imputed input data.
+        """
+        check_is_fitted(self, 'initial_imputer_')
+
+        X, X_filled, mask_missing_values = self._initial_imputation(X)
+
+        # edge case: in case the user specifies 0 for n_imputations,
+        # then there is no need to do burn in and the result should be
+        # just the initial imputation (before clipping)
+        if self.n_imputations < 1:
+            return X_filled
+
+        X_filled = np.clip(X_filled, self._min_value, self._max_value)
+
+        n_rounds = self.n_burn_in + self.n_imputations
+        n_imputations = len(self.imputation_sequence_)
+        imputations_per_round = n_imputations // n_rounds
+        i_rnd = 0
+        Xt = np.zeros(X.shape, dtype=X.dtype)
+        if self.verbose > 0:
+            print("[MICE] Completing matrix with shape %s" % (X.shape,))
+        start_t = time()
+        for it, predictor_triplet in enumerate(self.imputation_sequence_):
+            X_filled, _ = self._impute_one_feature(
+                X_filled,
+                mask_missing_values,
+                predictor_triplet.feat_idx,
+                predictor_triplet.neighbor_feat_idx,
+                predictor=predictor_triplet.predictor,
+                fit_mode=False
+            )
+            if not (it + 1) % imputations_per_round:
+                if i_rnd >= self.n_burn_in:
+                    Xt += X_filled
+                if self.verbose > 1:
+                    print('[MICE] Ending imputation round '
+                          '%d/%d, elapsed time %0.2f'
+                          % (i_rnd + 1, n_rounds, time() - start_t))
+                i_rnd += 1
+
+        Xt /= self.n_imputations
+        Xt[~mask_missing_values] = X[~mask_missing_values]
+        return Xt
+
+    def fit(self, X, y=None):
+        """Fits the imputer on X and return self.
+
+        Parameters
+        ----------
+        X : array-like, shape (n_samples, n_features)
+            Input data, where "n_samples" is the number of samples and
+            "n_features" is the number of features.
+
+        y : ignored
+
+        Returns
+        -------
+        self : object
+            Returns self.
+        """
+        self.fit_transform(X)
+        return self
diff --git a/sklearn/utils/estimator_checks.py b/sklearn/utils/estimator_checks.py
index 9a321e914b238..8508e166cd8f9 100644
--- a/sklearn/utils/estimator_checks.py
+++ b/sklearn/utils/estimator_checks.py
@@ -72,7 +72,7 @@
                 'OrthogonalMatchingPursuit', 'PLSCanonical', 'PLSRegression',
                 'RANSACRegressor', 'RadiusNeighborsRegressor',
                 'RandomForestRegressor', 'Ridge', 'RidgeCV']
-
+ALLOW_NAN = ['Imputer', 'SimpleImputer', 'MICEImputer']
 
 def _yield_non_meta_checks(name, estimator):
     yield check_estimators_dtypes
@@ -93,7 +93,7 @@ def _yield_non_meta_checks(name, estimator):
         # cross-decomposition's "transform" returns X and Y
         yield check_pipeline_consistency
 
-    if name not in ['SimpleImputer', 'Imputer']:
+    if name not in ALLOW_NAN:
         # Test that all estimators check their input for NaN's and infs
         yield check_estimators_nan_inf
 


Test patch:
diff --git a/sklearn/tests/test_impute.py b/sklearn/tests/test_impute.py
index f2bf5912e2213..954a016a835bb 100644
--- a/sklearn/tests/test_impute.py
+++ b/sklearn/tests/test_impute.py
@@ -1,14 +1,19 @@
+from __future__ import division
+
+import pytest
 
 import numpy as np
 from scipy import sparse
 
-from sklearn.utils.testing import assert_equal
+from sklearn.utils.testing import assert_allclose
 from sklearn.utils.testing import assert_array_equal
 from sklearn.utils.testing import assert_array_almost_equal
 from sklearn.utils.testing import assert_raises
 from sklearn.utils.testing import assert_false
 
-from sklearn.impute import SimpleImputer
+from sklearn.impute import SimpleImputer, MICEImputer
+from sklearn.dummy import DummyRegressor
+from sklearn.linear_model import BayesianRidge, ARDRegression
 from sklearn.pipeline import Pipeline
 from sklearn.model_selection import GridSearchCV
 from sklearn import tree
@@ -61,10 +66,14 @@ def test_imputation_shape():
 
     for strategy in ['mean', 'median', 'most_frequent']:
         imputer = SimpleImputer(strategy=strategy)
-        X_imputed = imputer.fit_transform(X)
-        assert_equal(X_imputed.shape, (10, 2))
         X_imputed = imputer.fit_transform(sparse.csr_matrix(X))
-        assert_equal(X_imputed.shape, (10, 2))
+        assert X_imputed.shape == (10, 2)
+        X_imputed = imputer.fit_transform(X)
+        assert X_imputed.shape == (10, 2)
+
+        mice_imputer = MICEImputer(initial_strategy=strategy)
+        X_imputed = mice_imputer.fit_transform(X)
+        assert X_imputed.shape == (10, 2)
 
 
 def safe_median(arr, *args, **kwargs):
@@ -257,3 +266,227 @@ def test_imputation_copy():
     Xt = imputer.fit(X).transform(X)
     Xt.data[0] = -1
     assert_false(np.all(X.data == Xt.data))
+
+    # Note: If X is sparse and if missing_values=0, then a (dense) copy of X is
+    # made, even if copy=False.
+
+
+def test_mice_rank_one():
+    rng = np.random.RandomState(0)
+    d = 100
+    A = rng.rand(d, 1)
+    B = rng.rand(1, d)
+    X = np.dot(A, B)
+    nan_mask = rng.rand(d, d) < 0.5
+    X_missing = X.copy()
+    X_missing[nan_mask] = np.nan
+
+    imputer = MICEImputer(n_imputations=5,
+                          n_burn_in=5,
+                          verbose=True,
+                          random_state=rng)
+    X_filled = imputer.fit_transform(X_missing)
+    assert_allclose(X_filled, X, atol=0.001)
+
+
+@pytest.mark.parametrize(
+    "imputation_order",
+    ['random', 'roman', 'ascending', 'descending', 'arabic']
+)
+def test_mice_imputation_order(imputation_order):
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 10
+    X = sparse_random_matrix(n, d, density=0.10, random_state=rng).toarray()
+    X[:, 0] = 1  # this column should not be discarded by MICEImputer
+
+    imputer = MICEImputer(missing_values=0,
+                          n_imputations=1,
+                          n_burn_in=1,
+                          n_nearest_features=5,
+                          min_value=0,
+                          max_value=1,
+                          verbose=False,
+                          imputation_order=imputation_order,
+                          random_state=rng)
+    imputer.fit_transform(X)
+    ordered_idx = [i.feat_idx for i in imputer.imputation_sequence_]
+    if imputation_order == 'roman':
+        assert np.all(ordered_idx[:d-1] == np.arange(1, d))
+    elif imputation_order == 'arabic':
+        assert np.all(ordered_idx[:d-1] == np.arange(d-1, 0, -1))
+    elif imputation_order == 'random':
+        ordered_idx_round_1 = ordered_idx[:d-1]
+        ordered_idx_round_2 = ordered_idx[d-1:]
+        assert ordered_idx_round_1 != ordered_idx_round_2
+    elif 'ending' in imputation_order:
+        assert len(ordered_idx) == 2 * (d - 1)
+
+
+@pytest.mark.parametrize(
+    "predictor",
+    [DummyRegressor(), BayesianRidge(), ARDRegression()]
+)
+def test_mice_predictors(predictor):
+    rng = np.random.RandomState(0)
+
+    n = 100
+    d = 10
+    X = sparse_random_matrix(n, d, density=0.10, random_state=rng).toarray()
+
+    imputer = MICEImputer(missing_values=0,
+                          n_imputations=1,
+                          n_burn_in=1,
+                          predictor=predictor,
+                          random_state=rng)
+    imputer.fit_transform(X)
+
+    # check that types are correct for predictors
+    hashes = []
+    for triplet in imputer.imputation_sequence_:
+        assert triplet.predictor
+        hashes.append(id(triplet.predictor))
+
+    # check that each predictor is unique
+    assert len(set(hashes)) == len(hashes)
+
+
+def test_mice_clip():
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 10
+    X = sparse_random_matrix(n, d, density=0.10,
+                             random_state=rng).toarray()
+
+    imputer = MICEImputer(missing_values=0,
+                          n_imputations=1,
+                          n_burn_in=1,
+                          min_value=0.1,
+                          max_value=0.2,
+                          random_state=rng)
+
+    Xt = imputer.fit_transform(X)
+    assert_allclose(np.min(Xt[X == 0]), 0.1)
+    assert_allclose(np.max(Xt[X == 0]), 0.2)
+    assert_allclose(Xt[X != 0], X[X != 0])
+
+
+@pytest.mark.parametrize(
+    "strategy",
+    ["mean", "median", "most_frequent"]
+)
+def test_mice_missing_at_transform(strategy):
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 10
+    X_train = rng.randint(low=0, high=3, size=(n, d))
+    X_test = rng.randint(low=0, high=3, size=(n, d))
+
+    X_train[:, 0] = 1  # definitely no missing values in 0th column
+    X_test[0, 0] = 0  # definitely missing value in 0th column
+
+    mice = MICEImputer(missing_values=0,
+                       n_imputations=1,
+                       n_burn_in=1,
+                       initial_strategy=strategy,
+                       random_state=rng).fit(X_train)
+    initial_imputer = SimpleImputer(missing_values=0,
+                                    strategy=strategy).fit(X_train)
+
+    # if there were no missing values at time of fit, then mice will
+    # only use the initial imputer for that feature at transform
+    assert np.all(mice.transform(X_test)[:, 0] ==
+                  initial_imputer.transform(X_test)[:, 0])
+
+
+def test_mice_transform_stochasticity():
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 10
+    X = sparse_random_matrix(n, d, density=0.10,
+                             random_state=rng).toarray()
+
+    imputer = MICEImputer(missing_values=0,
+                          n_imputations=1,
+                          n_burn_in=1,
+                          random_state=rng)
+    imputer.fit(X)
+
+    X_fitted_1 = imputer.transform(X)
+    X_fitted_2 = imputer.transform(X)
+
+    # sufficient to assert that the means are not the same
+    assert np.mean(X_fitted_1) != pytest.approx(np.mean(X_fitted_2))
+
+
+def test_mice_no_missing():
+    rng = np.random.RandomState(0)
+    X = rng.rand(100, 100)
+    X[:, 0] = np.nan
+    m1 = MICEImputer(n_imputations=10, random_state=rng)
+    m2 = MICEImputer(n_imputations=10, random_state=rng)
+    pred1 = m1.fit(X).transform(X)
+    pred2 = m2.fit_transform(X)
+    # should exclude the first column entirely
+    assert_allclose(X[:, 1:], pred1)
+    # fit and fit_transform should both be identical
+    assert_allclose(pred1, pred2)
+
+
+@pytest.mark.parametrize(
+    "rank",
+    [3, 5]
+)
+def test_mice_transform_recovery(rank):
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 100
+    A = rng.rand(n, rank)
+    B = rng.rand(rank, d)
+    X_filled = np.dot(A, B)
+    # half is randomly missing
+    nan_mask = rng.rand(n, d) < 0.5
+    X_missing = X_filled.copy()
+    X_missing[nan_mask] = np.nan
+
+    # split up data in half
+    n = n // 2
+    X_train = X_missing[:n]
+    X_test_filled = X_filled[n:]
+    X_test = X_missing[n:]
+
+    imputer = MICEImputer(n_imputations=10,
+                          n_burn_in=10,
+                          verbose=True,
+                          random_state=rng).fit(X_train)
+    X_test_est = imputer.transform(X_test)
+    assert_allclose(X_test_filled, X_test_est, rtol=1e-5, atol=0.1)
+
+
+def test_mice_additive_matrix():
+    rng = np.random.RandomState(0)
+    n = 100
+    d = 10
+    A = rng.randn(n, d)
+    B = rng.randn(n, d)
+    X_filled = np.zeros(A.shape)
+    for i in range(d):
+        for j in range(d):
+            X_filled[:, (i+j) % d] += (A[:, i] + B[:, j]) / 2
+    # a quarter is randomly missing
+    nan_mask = rng.rand(n, d) < 0.25
+    X_missing = X_filled.copy()
+    X_missing[nan_mask] = np.nan
+
+    # split up data
+    n = n // 2
+    X_train = X_missing[:n]
+    X_test_filled = X_filled[n:]
+    X_test = X_missing[n:]
+
+    imputer = MICEImputer(n_imputations=25,
+                          n_burn_in=10,
+                          verbose=True,
+                          random_state=rng).fit(X_train)
+    X_test_est = imputer.transform(X_test)
+    assert_allclose(X_test_filled, X_test_est, atol=0.01)