diff --git a/RELEASE-NOTES.md b/RELEASE-NOTES.md
index 50f95be271..6fc30bad59 100644
--- a/RELEASE-NOTES.md
+++ b/RELEASE-NOTES.md
@@ -18,6 +18,7 @@
 - Deprecated `sd` in version 3.7 has been replaced by `sigma` now raises `DepreciationWarning` on using `sd` in continuous, mixed and timeseries distributions. (see #3837 and #3688).
 - In named models, `pm.Data` objects now get model-relative names (see [#3843](https://github.com/pymc-devs/pymc3/pull/3843)).
 - `pm.sample` now takes 1000 draws and 1000 tuning samples by default, instead of 500 previously (see [#3855](https://github.com/pymc-devs/pymc3/pull/3855)).
+- Dropped the outdated 'nuts' initialization method for `pm.sample` (see [#3863](https://github.com/pymc-devs/pymc3/pull/3863)).
 
 ## PyMC3 3.8 (November 29 2019)
 
diff --git a/pymc3/sampling.py b/pymc3/sampling.py
index 8ecce6828b..eab382ce87 100644
--- a/pymc3/sampling.py
+++ b/pymc3/sampling.py
@@ -274,7 +274,6 @@ def sample(
         * advi: Run ADVI to estimate posterior mean and diagonal mass matrix.
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map: Use the MAP as starting point. This is discouraged.
-        * nuts: Run NUTS and estimate posterior mean and mass matrix from the trace.
         * adapt_full: Adapt a dense mass matrix using the sample covariances
     step: function or iterable of functions
         A step function or collection of functions. If there are variables without step methods,
@@ -282,8 +281,7 @@ def sample(
         method will be used, if appropriate to the model; this is a good default for beginning
         users.
     n_init: int
-        Number of iterations of initializer. Only works for 'nuts' and 'ADVI'.
-        If 'ADVI', number of iterations, if 'nuts', number of draws.
+        Number of iterations of initializer. Only works for 'ADVI' init methods.
     start: dict, or array of dict
         Starting point in parameter space (or partial point)
         Defaults to ``trace.point(-1))`` if there is a trace provided and model.test_point if not
@@ -1865,14 +1863,11 @@ def init_nuts(
         * advi: Run ADVI to estimate posterior mean and diagonal mass matrix.
         * advi_map: Initialize ADVI with MAP and use MAP as starting point.
         * map: Use the MAP as starting point. This is discouraged.
-        * nuts: Run NUTS and estimate posterior mean and mass matrix from
-          the trace.
         * adapt_full: Adapt a dense mass matrix using the sample covariances
     chains: int
         Number of jobs to start.
     n_init: int
-        Number of iterations of initializer
-        If 'ADVI', number of iterations, if 'nuts', number of draws.
+        Number of iterations of initializer. Only works for 'ADVI' init methods.
     model: Model (optional if in ``with`` context)
     progressbar: bool
         Whether or not to display a progressbar for advi sampling.
@@ -2001,13 +1996,6 @@ def init_nuts(
         cov = pm.find_hessian(point=start)
         start = [start] * chains
         potential = quadpotential.QuadPotentialFull(cov)
-    elif init == "nuts":
-        init_trace = pm.sample(
-            draws=n_init, step=pm.NUTS(), tune=n_init // 2, random_seed=random_seed
-        )
-        cov = np.atleast_1d(pm.trace_cov(init_trace))
-        start = list(np.random.choice(init_trace, chains))
-        potential = quadpotential.QuadPotentialFull(cov)
     elif init == "adapt_full":
         start = [model.test_point] * chains
         mean = np.mean([model.dict_to_array(vals) for vals in start], axis=0)
diff --git a/pymc3/tests/test_sampling.py b/pymc3/tests/test_sampling.py
index 6c0eece037..527c944ab7 100644
--- a/pymc3/tests/test_sampling.py
+++ b/pymc3/tests/test_sampling.py
@@ -87,7 +87,7 @@ def test_sample(self):
 
     def test_sample_init(self):
         with self.model:
-            for init in ("advi", "advi_map", "map", "nuts"):
+            for init in ("advi", "advi_map", "map"):
                 pm.sample(
                     init=init,
                     tune=0,
@@ -675,7 +675,6 @@ def test_sample_posterior_predictive_w(self):
         "advi+adapt_diag_grad",
         "map",
         "advi_map",
-        "nuts",
     ],
 )
 def test_exec_nuts_init(method):