diff --git a/doc/whats_new/v1.0.rst b/doc/whats_new/v1.0.rst
index 382ff363e0db7..004762498d49e 100644
--- a/doc/whats_new/v1.0.rst
+++ b/doc/whats_new/v1.0.rst
@@ -115,6 +115,14 @@ Changelog
   class methods and will be removed in 1.2.
   :pr:`18543` by `Guillaume Lemaitre`_.
 
+:mod:`sklearn.multioutput`
+..........................
+
+- |Fix| :func:`multioutput.ClassifierChain.decision_function` and
+  :func:`multioutput.ClassifierChain.predict_proba` now both return a list of
+  ``n_outputs`` arrays of shape `(n_samples, n_classes).` :pr:`14654` by
+  :user:`Agamemnon Krasoulis <agamemnonc>`.
+
 :mod:`sklearn.naive_bayes`
 ..........................
 
diff --git a/sklearn/multioutput.py b/sklearn/multioutput.py
index 9987c01b13187..27c143074cd72 100644
--- a/sklearn/multioutput.py
+++ b/sklearn/multioutput.py
@@ -638,9 +638,15 @@ class labels for each estimator in the chain.
            [1., 0., 0.],
            [0., 1., 0.]])
     >>> chain.predict_proba(X_test)
-    array([[0.8387..., 0.9431..., 0.4576...],
-           [0.8878..., 0.3684..., 0.2640...],
-           [0.0321..., 0.9935..., 0.0625...]])
+    [array([[0.16126878, 0.83873122],
+        [0.11218344, 0.88781656],
+        [0.96786386, 0.03213614]]),
+    array([[0.05685769, 0.94314231],
+        [0.6315953 , 0.3684047 ],
+        [0.00640331, 0.99359669]]),
+    array([[0.5423851 , 0.4576149 ],
+        [0.73590132, 0.26409868],
+        [0.93742079, 0.06257921]])]
 
     See Also
     --------
@@ -684,10 +690,13 @@ def predict_proba(self, X):
 
         Returns
         -------
-        Y_prob : array-like of shape (n_samples, n_classes)
+        Y_prob : list of n_outputs ndarray of shape (n_samples, n_classes)
+            The class probabilities of the input samples. The order of the
+            classes for each output corresponds to the respective entry of
+            the attribute `classes_`.
         """
         X = check_array(X, accept_sparse=True)
-        Y_prob_chain = np.zeros((X.shape[0], len(self.estimators_)))
+        Y_prob_chain = []
         Y_pred_chain = np.zeros((X.shape[0], len(self.estimators_)))
         for chain_idx, estimator in enumerate(self.estimators_):
             previous_predictions = Y_pred_chain[:, :chain_idx]
@@ -695,11 +704,11 @@ def predict_proba(self, X):
                 X_aug = sp.hstack((X, previous_predictions))
             else:
                 X_aug = np.hstack((X, previous_predictions))
-            Y_prob_chain[:, chain_idx] = estimator.predict_proba(X_aug)[:, 1]
+            Y_prob_chain.append(estimator.predict_proba(X_aug))
             Y_pred_chain[:, chain_idx] = estimator.predict(X_aug)
         inv_order = np.empty_like(self.order_)
         inv_order[self.order_] = np.arange(len(self.order_))
-        Y_prob = Y_prob_chain[:, inv_order]
+        Y_prob = [Y_prob_chain[i] for i in inv_order]
 
         return Y_prob
 
@@ -713,11 +722,12 @@ def decision_function(self, X):
 
         Returns
         -------
-        Y_decision : array-like of shape (n_samples, n_classes)
-            Returns the decision function of the sample for each model
-            in the chain.
+        Y_decision : list of n_outputs ndarray of shape (n_samples, n_classes)
+            Decision function of the input samples for each model
+            in the chain.  The order of the classes for each output corresponds
+            to the respective entry of the attribute `classes_`.
         """
-        Y_decision_chain = np.zeros((X.shape[0], len(self.estimators_)))
+        Y_decision_chain = []
         Y_pred_chain = np.zeros((X.shape[0], len(self.estimators_)))
         for chain_idx, estimator in enumerate(self.estimators_):
             previous_predictions = Y_pred_chain[:, :chain_idx]
@@ -725,18 +735,17 @@ def decision_function(self, X):
                 X_aug = sp.hstack((X, previous_predictions))
             else:
                 X_aug = np.hstack((X, previous_predictions))
-            Y_decision_chain[:, chain_idx] = estimator.decision_function(X_aug)
+            Y_decision_chain.append(estimator.decision_function(X_aug))
             Y_pred_chain[:, chain_idx] = estimator.predict(X_aug)
 
         inv_order = np.empty_like(self.order_)
         inv_order[self.order_] = np.arange(len(self.order_))
-        Y_decision = Y_decision_chain[:, inv_order]
+        Y_decision = [Y_decision_chain[i] for i in inv_order]
 
         return Y_decision
 
     def _more_tags(self):
-        return {'_skip_test': True,
-                'multioutput_only': True}
+        return {'_skip_test': True}
 
 
 class RegressorChain(MetaEstimatorMixin, RegressorMixin, _BaseChain):
diff --git a/sklearn/tests/test_multioutput.py b/sklearn/tests/test_multioutput.py
index edfcdef1bf89c..65fabd0cd1b7a 100644
--- a/sklearn/tests/test_multioutput.py
+++ b/sklearn/tests/test_multioutput.py
@@ -424,8 +424,9 @@ def test_classifier_chain_fit_and_predict_with_linear_svc():
     assert Y_pred.shape == Y.shape
 
     Y_decision = classifier_chain.decision_function(X)
+    Y_binary = [Y_decision[i] >= 0 for i in range(Y.shape[1])]
+    Y_binary = np.asarray(Y_binary).T
 
-    Y_binary = (Y_decision >= 0)
     assert_array_equal(Y_binary, Y_pred)
     assert not hasattr(classifier_chain, 'predict_proba')
 
@@ -481,7 +482,8 @@ def test_base_chain_fit_and_predict():
                 list(range(X.shape[1], X.shape[1] + Y.shape[1])))
 
     Y_prob = chains[1].predict_proba(X)
-    Y_binary = (Y_prob >= .5)
+    Y_binary = [np.argmax(Y_prob[i], axis=1) for i in range(Y.shape[1])]
+    Y_binary = np.asarray(Y_binary).T
     assert_array_equal(Y_binary, Y_pred)
 
     assert isinstance(chains[1], ClassifierMixin)