Repository: pymc-devs/pymc

Issue:
Interpolated distribution in V4
## Description of your problem

While working on the issue https://github.com/pymc-devs/pymc-examples/issues/115, where I am updating the code from PyMC v3 to v4, I got this error:

```python
File ~/Documents/tech_world/bayesian/pymc/pymc/model.py:1790, in Model.check_start_vals(self, start)
   1787 initial_eval = self.point_logps(point=elem)
   1789 if not all(np.isfinite(v) for v in initial_eval.values()):
-> 1790     raise SamplingError(
   1791         "Initial evaluation of model at starting point failed!\n"
   1792         f"Starting values:\n{elem}\n\n"
...
Starting values:
{'alpha_interval__': array(nan), 'beta0_interval__': array(nan), 'beta1_interval__': array(nan)}

Initial evaluation results:
{'alpha': nan, 'beta0': nan, 'beta1': nan, 'Y_obs': nan}
```
Trying to find the reason, I realized that the next basic model works in **v3 but not in V4.**

**Please provide a minimal, self-contained, and reproducible example.**
```python
import numpy as np
import pymc as pm

from scipy import stats

model = pm.Model()

x_points = np.linspace(0, 10, 100)
pdf_points = stats.norm.pdf(x_points, loc=1, scale=1)
with model:    
    # alpha = from_posterior("alpha", trace["alpha"])
    alpha = pm.Interpolated("alpha", x_points, pdf_points)
    # Expected value of outcome
    pm.sample(1000)
```

**Please provide the full traceback.**

<details><summary>Complete error traceback</summary>

```python
---------------------------------------------------------------------------
SamplingError                             Traceback (most recent call last)
/Users/symeneses/Documents/tech_world/bayesian/pymc-examples/examples/howto/updating_priors.ipynb Cell 16' in <cell line: 10>()
     [12](vscode-notebook-cell:/Users/symeneses/Documents/tech_world/bayesian/pymc-examples/examples/howto/updating_priors.ipynb#ch0000023?line=11) alpha = pm.Interpolated("alpha", x_points, pdf_points)
     [13](vscode-notebook-cell:/Users/symeneses/Documents/tech_world/bayesian/pymc-examples/examples/howto/updating_priors.ipynb#ch0000023?line=12) # Expected value of outcome
---> [14](vscode-notebook-cell:/Users/symeneses/Documents/tech_world/bayesian/pymc-examples/examples/howto/updating_priors.ipynb#ch0000023?line=13) pm.sample(1000)

File ~/Documents/tech_world/bayesian/pymc/pymc/sampling.py:558, in sample(draws, step, init, n_init, initvals, trace, chain_idx, chains, cores, tune, progressbar, model, random_seed, discard_tuned_samples, compute_convergence_checks, callback, jitter_max_retries, return_inferencedata, idata_kwargs, mp_ctx, **kwargs)
    556 # One final check that shapes and logps at the starting points are okay.
    557 for ip in initial_points:
--> 558     model.check_start_vals(ip)
    559     _check_start_shape(model, ip)
    561 sample_args = {
    562     "draws": draws,
    563     "step": step,
   (...)
    573     "discard_tuned_samples": discard_tuned_samples,
    574 }

File ~/Documents/tech_world/bayesian/pymc/pymc/model.py:1790, in Model.check_start_vals(self, start)
   1787 initial_eval = self.point_logps(point=elem)
   1789 if not all(np.isfinite(v) for v in initial_eval.values()):
-> 1790     raise SamplingError(
   1791         "Initial evaluation of model at starting point failed!\n"
   1792         f"Starting values:\n{elem}\n\n"
```

</details>

**Please provide any additional information below.**

Changing the number of sample in `linspace` from 100 to 10 , it works also in V4.
```python
import numpy as np
import pymc as pm

from scipy import stats

model = pm.Model()

x_points = np.linspace(0, 10, 10)
pdf_points = stats.norm.pdf(x_points, loc=1, scale=1)
with model:    
    # alpha = from_posterior("alpha", trace["alpha"])
    alpha = pm.Interpolated("alpha", x_points, pdf_points)
    # Expected value of outcome
    pm.sample(1000)
```

## Versions and main components

* PyMC/PyMC3 Version: v4.0.1 (It works in 3.11.4)
* Aesara/Theano Version: 2.7.3
* Python Version: 3.10
* Operating system: macOS 12.4
* How did you install PyMC/PyMC3: conda


Discussion:
codecov[bot]: # [Codecov](https://codecov.io/gh/pymc-devs/pymc/pull/5986?src=pr&el=h1&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) Report
> Merging [#5986](https://codecov.io/gh/pymc-devs/pymc/pull/5986?src=pr&el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) (bca54ad) into [main](https://codecov.io/gh/pymc-devs/pymc/commit/2a9e86c641ccd2d8da87268eefc3bbb72d21fbe5?el=desc&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) (2a9e86c) will **decrease** coverage by `14.96%`.
> The diff coverage is `0.00%`.

[![Impacted file tree graph](https://codecov.io/gh/pymc-devs/pymc/pull/5986/graphs/tree.svg?width=650&height=150&src=pr&token=JFuXtOJ4Cb&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs)](https://codecov.io/gh/pymc-devs/pymc/pull/5986?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs)

```diff
@@             Coverage Diff             @@
##             main    #5986       +/-   ##
===========================================
- Coverage   87.31%   72.34%   -14.97%     
===========================================
  Files          73       73               
  Lines       13253    13254        +1     
===========================================
- Hits        11572     9589     -1983     
- Misses       1681     3665     +1984     
```


| [Impacted Files](https://codecov.io/gh/pymc-devs/pymc/pull/5986?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) | Coverage Δ | |
|---|---|---|
| [pymc/distributions/continuous.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL2NvbnRpbnVvdXMucHk=) | `74.52% <0.00%> (-23.27%)` | :arrow_down: |
| [pymc/smc/smc.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9zbWMvc21jLnB5) | `18.89% <0.00%> (-77.56%)` | :arrow_down: |
| [pymc/tuning/starting.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy90dW5pbmcvc3RhcnRpbmcucHk=) | `19.32% <0.00%> (-73.11%)` | :arrow_down: |
| [pymc/variational/inference.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy92YXJpYXRpb25hbC9pbmZlcmVuY2UucHk=) | `19.25% <0.00%> (-66.85%)` | :arrow_down: |
| [pymc/smc/sample\_smc.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9zbWMvc2FtcGxlX3NtYy5weQ==) | `16.77% <0.00%> (-66.45%)` | :arrow_down: |
| [pymc/func\_utils.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9mdW5jX3V0aWxzLnB5) | `22.50% <0.00%> (-65.00%)` | :arrow_down: |
| [pymc/variational/test\_functions.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy92YXJpYXRpb25hbC90ZXN0X2Z1bmN0aW9ucy5weQ==) | `38.09% <0.00%> (-61.91%)` | :arrow_down: |
| [pymc/distributions/simulator.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL3NpbXVsYXRvci5weQ==) | `31.25% <0.00%> (-56.25%)` | :arrow_down: |
| [pymc/distributions/bound.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy9kaXN0cmlidXRpb25zL2JvdW5kLnB5) | `45.54% <0.00%> (-54.46%)` | :arrow_down: |
| [pymc/variational/stein.py](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs#diff-cHltYy92YXJpYXRpb25hbC9zdGVpbi5weQ==) | `48.33% <0.00%> (-51.67%)` | :arrow_down: |
| ... and [34 more](https://codecov.io/gh/pymc-devs/pymc/pull/5986/diff?src=pr&el=tree-more&utm_medium=referral&utm_source=github&utm_content=comment&utm_campaign=pr+comments&utm_term=pymc-devs) | |


Patch:
diff --git a/.pre-commit-config.yaml b/.pre-commit-config.yaml
index 2db19a5b4..e0118cd33 100644
--- a/.pre-commit-config.yaml
+++ b/.pre-commit-config.yaml
@@ -26,7 +26,7 @@ repos:
       - types-filelock
       - types-setuptools
       - arviz
-      - aesara==2.7.7
+      - aesara==2.7.8
       - aeppl==0.0.32
       always_run: true
       require_serial: true
diff --git a/conda-envs/environment-dev.yml b/conda-envs/environment-dev.yml
index 4bd145aac..aa0e56214 100644
--- a/conda-envs/environment-dev.yml
+++ b/conda-envs/environment-dev.yml
@@ -6,7 +6,7 @@ channels:
 dependencies:
 # Base dependencies
 - aeppl=0.0.32
-- aesara=2.7.7
+- aesara=2.7.8
 - arviz>=0.12.0
 - blas
 - cachetools>=4.2.1
diff --git a/conda-envs/windows-environment-dev.yml b/conda-envs/windows-environment-dev.yml
index 198321e47..ccfe32206 100644
--- a/conda-envs/windows-environment-dev.yml
+++ b/conda-envs/windows-environment-dev.yml
@@ -6,7 +6,7 @@ channels:
 dependencies:
 # Base dependencies (see install guide for Windows)
 - aeppl=0.0.32
-- aesara=2.7.7
+- aesara=2.7.8
 - arviz>=0.12.0
 - blas
 - cachetools>=4.2.1
diff --git a/docs/source/api/samplers.rst b/docs/source/api/samplers.rst
index 1cd43b1da..0acbc9df2 100644
--- a/docs/source/api/samplers.rst
+++ b/docs/source/api/samplers.rst
@@ -52,17 +52,6 @@ Metropolis family
     PoissonProposal
     UniformProposal
 
-MLDA family
------------
-
-.. autosummary::
-   :toctree: generated/
-
-   MLDA
-   DEMetropolisZMLDA
-   MetropolisMLDA
-   RecursiveDAProposal
-
 Other step methods
 ------------------
 
diff --git a/docs/source/contributing/developer_guide.md b/docs/source/contributing/developer_guide.md
index 5137e1fd3..f40c14621 100644
--- a/docs/source/contributing/developer_guide.md
+++ b/docs/source/contributing/developer_guide.md
@@ -151,7 +151,7 @@ def logp(self, value):
 ```
 
 In the ``logp`` method, parameters and values are either Aesara tensors, or could be converted to tensors.
-It is rather convenient as the evaluation of logp is represented as a tensor (``RV.logpt``), and when we linked different ``logp`` together (e.g., summing all ``RVs.logpt`` to get the model totall logp) the dependence  is taken care of by Aesara when the graph is built and compiled.
+It is rather convenient as the evaluation of logp is represented as a tensor (``RV.logpt``), and when we linked different ``logp`` together (e.g., summing all ``RVs.logpt`` to get the model total logp) the dependence  is taken care of by Aesara when the graph is built and compiled.
 Again, since the compiled function depends on the nodes that already in the graph, whenever you want to generate a new function that takes new input tensors you either need to regenerate the graph with the appropriate dependencies, or replace the node by editing the existing graph.
 In PyMC we use the second approach by using ``aesara.clone_replace()`` when it is needed.
 
diff --git a/docs/source/contributing/jupyter_style.md b/docs/source/contributing/jupyter_style.md
index d4ef53584..de34ba8db 100644
--- a/docs/source/contributing/jupyter_style.md
+++ b/docs/source/contributing/jupyter_style.md
@@ -53,7 +53,7 @@ This guide does not teach nor cover MyST extensively, only gives some opinionate
 * Above all, stay consistent with variable names within the notebook. Notebooks using multiple names for the same variable will not be merged.
 * Use meaningful variable names wherever possible. Our users come from different backgrounds and not everyone is familiar with the same naming conventions.
   - Annotate dimensions too. Notebooks are published to be read, so even if the shape is derived
-    from the inputs or you don't like to use named dims and don't use them in your presonal
+    from the inputs or you don't like to use named dims and don't use them in your personal
     code, notebooks must use dims, even if annotating and not setting the shape.
     It makes the code easier to follow, especially for newcomers.
 * Sometimes it makes sense to use Greek letters to refer to variables, for example when writing equations, as this makes them easier to read. In that case, use LaTeX to insert the Greek letter like this `$\theta$` instead of using Unicode like `θ`.
@@ -397,6 +397,8 @@ After the notebook content finishes, there should be an `## Authors` section wit
 to provide attribution to the people who contributed to the the general pattern should be:
 
 ```markdown
+## Authors
+
 * <verb> by <author> on <date> ([repo#PR](https://link-to.pr))
 ```
 
@@ -413,6 +415,8 @@ adapted
   in order to include a link or reference to the original source:
 
   ```markdown
+  ## Authors
+
   Adapted from Alice's [blogpost](blog.alice.com) by Bob and Carol on ...
   ```
 
@@ -427,6 +431,8 @@ updated
 some examples:
 
 ```markdown
+## Authors
+
 * Authored by Chris Fonnesbeck in May, 2017 ([pymc#2124](https://github.com/pymc-devs/pymc/pull/2124))
 * Updated by Colin Carroll in June, 2018 ([pymc#3049](https://github.com/pymc-devs/pymc/pull/3049))
 * Updated by Alex Andorra in January, 2020 ([pymc#3765](https://github.com/pymc-devs/pymc/pull/3765))
@@ -437,6 +443,8 @@ some examples:
 and
 
 ```markdown
+## Authors
+
 * Adapted from chapter 5 of Bayesian Data Analysis 3rd Edition {cite:p}`gelman2013bayesian`
   by Demetri Pananos and Junpeng Lao on June, 2018 ([pymc#3054](https://github.com/pymc-devs/pymc/pull/3054))
 * Reexecuted by Ravin Kumar with PyMC 3.6 on March, 2019 ([pymc#3397](https://github.com/pymc-devs/pymc/pull/3397))
diff --git a/pymc/distributions/continuous.py b/pymc/distributions/continuous.py
index 2de42a399..0dce86078 100644
--- a/pymc/distributions/continuous.py
+++ b/pymc/distributions/continuous.py
@@ -3734,8 +3734,11 @@ class Interpolated(BoundedContinuous):
         return super().dist([x_points, pdf_points, cdf_points], **kwargs)
 
     def moment(rv, size, x_points, pdf_points, cdf_points):
-        # cdf_points argument is unused
-        moment = at.sum(at.mul(x_points, pdf_points))
+        """
+        Estimates the expectation integral using the trapezoid rule; cdf_points are not used.
+        """
+        x_fx = at.mul(x_points, pdf_points)  # x_i * f(x_i) for all xi's in x_points
+        moment = at.sum(at.mul(at.diff(x_points), x_fx[1:] + x_fx[:-1])) / 2
 
         if not rv_size_is_none(size):
             moment = at.full(size, moment)
diff --git a/pymc/step_methods/__init__.py b/pymc/step_methods/__init__.py
index 8eb79a44b..d7beea111 100644
--- a/pymc/step_methods/__init__.py
+++ b/pymc/step_methods/__init__.py
@@ -28,12 +28,6 @@ from pymc.step_methods.metropolis import (
     PoissonProposal,
     UniformProposal,
 )
-from pymc.step_methods.mlda import (
-    MLDA,
-    DEMetropolisZMLDA,
-    MetropolisMLDA,
-    RecursiveDAProposal,
-)
 from pymc.step_methods.slicer import Slice
 
 STEP_METHODS = (
diff --git a/pymc/step_methods/hmc/hmc.py b/pymc/step_methods/hmc/hmc.py
index 096f0ea29..6f1522596 100644
--- a/pymc/step_methods/hmc/hmc.py
+++ b/pymc/step_methods/hmc/hmc.py
@@ -77,7 +77,7 @@ class HamiltonianMC(BaseHMC):
         scaling: array_like, ndim = {1,2}
             The inverse mass, or precision matrix. One dimensional arrays are
             interpreted as diagonal matrices. If `is_cov` is set to True,
-            this will be interpreded as the mass or covariance matrix.
+            this will be interpreted as the mass or covariance matrix.
         is_cov: bool, default=False
             Treat the scaling as mass or covariance matrix.
         potential: Potential, optional
diff --git a/pymc/step_methods/hmc/nuts.py b/pymc/step_methods/hmc/nuts.py
index 49906f6bc..c9005c964 100644
--- a/pymc/step_methods/hmc/nuts.py
+++ b/pymc/step_methods/hmc/nuts.py
@@ -158,7 +158,7 @@ class NUTS(BaseHMC):
         scaling: array_like, ndim = {1,2}
             The inverse mass, or precision matrix. One dimensional arrays are
             interpreted as diagonal matrices. If `is_cov` is set to True,
-            this will be interpreded as the mass or covariance matrix.
+            this will be interpreted as the mass or covariance matrix.
         is_cov: bool, default=False
             Treat the scaling as mass or covariance matrix.
         potential: Potential, optional
diff --git a/pymc/step_methods/mlda.py b/pymc/step_methods/mlda.py
deleted file mode 100644
index 5b6f35073..000000000
--- a/pymc/step_methods/mlda.py
+++ /dev/null
@@ -1,1099 +0,0 @@
-#   Copyright 2020 The PyMC Developers
-#
-#   Licensed under the Apache License, Version 2.0 (the "License");
-#   you may not use this file except in compliance with the License.
-#   You may obtain a copy of the License at
-#
-#       http://www.apache.org/licenses/LICENSE-2.0
-#
-#   Unless required by applicable law or agreed to in writing, software
-#   distributed under the License is distributed on an "AS IS" BASIS,
-#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
-#   See the License for the specific language governing permissions and
-#   limitations under the License.
-
-import logging
-import warnings
-
-from typing import Any, Dict, List, Optional, Tuple, Type, Union
-
-import arviz as az
-import numpy as np
-
-from aesara.tensor.sharedvar import TensorSharedVariable
-
-import pymc as pm
-
-from pymc.blocking import DictToArrayBijection, RaveledVars
-from pymc.model import Model, Point
-from pymc.step_methods.arraystep import ArrayStepShared, Competence, metrop_select
-from pymc.step_methods.compound import CompoundStep
-from pymc.step_methods.metropolis import DEMetropolisZ, Metropolis, Proposal, delta_logp
-
-__all__ = [
-    "MetropolisMLDA",
-    "DEMetropolisZMLDA",
-    "RecursiveDAProposal",
-    "MLDA",
-    "extract_Q_estimate",
-]
-
-
-class MetropolisMLDA(Metropolis):
-    """
-    Metropolis-Hastings sampling step tailored for use as base sampler in MLDA.
-    """
-
-    name = "metropolis_mlda"
-
-    def __init__(self, *args, **kwargs):
-        """
-        Initialise MetropolisMLDA. This is a mix of the parent's class' initialisation
-        and some extra code specific for MLDA.
-        """
-        model = pm.modelcontext(kwargs.get("model", None))
-        initial_values = model.initial_point()
-
-        # flag to that variance reduction is activated - forces MetropolisMLDA
-        # to store quantities of interest in a register if True
-        self.mlda_variance_reduction = kwargs.pop("mlda_variance_reduction", False)
-        if self.mlda_variance_reduction:
-            # Subsampling rate of MLDA sampler one level up
-            self.mlda_subsampling_rate_above = kwargs.pop("mlda_subsampling_rate_above")
-            self.sub_counter = 0
-            self.Q_last = np.nan
-            self.Q_reg = [np.nan] * self.mlda_subsampling_rate_above
-
-        # call parent class __init__
-        super().__init__(*args, **kwargs)
-
-        # modify the delta function and point to model if VR is used
-        if self.mlda_variance_reduction:
-            self.model = model
-            self.delta_logp_factory = self.delta_logp
-            self.delta_logp = lambda q, q0: -self.delta_logp_factory(q0, q)
-
-    def reset_tuning(self):
-        """
-        Does not reset sampler parameters. Allows continuation with
-        the same settings when MetropolisMLDA steps are done in chunks
-        under MLDA.
-        """
-        return
-
-    def astep(self, q0):
-
-        q_new, stats = super().astep(q0)
-
-        # Add variance reduction functionality.
-        if self.mlda_variance_reduction:
-            if stats[0]["accepted"]:
-                self.Q_last = self.model.Q.get_value()
-            if self.sub_counter == self.mlda_subsampling_rate_above:
-                self.sub_counter = 0
-            self.Q_reg[self.sub_counter] = self.Q_last
-            self.sub_counter += 1
-
-        return q_new, stats
-
-
-class DEMetropolisZMLDA(DEMetropolisZ):
-    """
-    DEMetropolisZ sampling step tailored for use as base sampler in MLDA
-    """
-
-    name = "DEMetropolisZ_mlda"
-
-    def __init__(self, *args, **kwargs):
-        """
-        Initialise DEMetropolisZMLDA, uses parent class __init__
-        and extra code specific for use within MLDA.
-        """
-
-        # flag used for signaling the end of tuning
-        self.tuning_end_trigger = False
-
-        model = pm.modelcontext(kwargs.get("model", None))
-        initial_values = model.initial_point()
-
-        # flag to that variance reduction is activated - forces DEMetropolisZMLDA
-        # to store quantities of interest in a register if True
-        self.mlda_variance_reduction = kwargs.pop("mlda_variance_reduction", False)
-        if self.mlda_variance_reduction:
-            # Subsampling rate of MLDA sampler one level up
-            self.mlda_subsampling_rate_above = kwargs.pop("mlda_subsampling_rate_above")
-            self.sub_counter = 0
-            self.Q_last = np.nan
-            self.Q_reg = [np.nan] * self.mlda_subsampling_rate_above
-
-        # call parent class __init__
-        super().__init__(*args, **kwargs)
-
-        # modify the delta function and point to model if VR is used
-        if self.mlda_variance_reduction:
-            self.model = model
-            self.delta_logp_factory = self.delta_logp
-            self.delta_logp = lambda q, q0: -self.delta_logp_factory(q0, q)
-
-    def reset_tuning(self):
-        """Skips resetting of tuned sampler parameters
-        and history to their initial values. Allows
-        continuation with the same settings when
-        DEMetropolisZMLDA steps are done in chunks
-        under MLDA."""
-        return
-
-    def astep(self, q0):
-
-        q_new, stats = super().astep(q0)
-
-        # Add variance reduction functionality.
-        if self.mlda_variance_reduction:
-            if stats[0]["accepted"]:
-                self.Q_last = self.model.Q.get_value()
-            if self.sub_counter == self.mlda_subsampling_rate_above:
-                self.sub_counter = 0
-            self.Q_reg[self.sub_counter] = self.Q_last
-            self.sub_counter += 1
-
-        return q_new, stats
-
-    def stop_tuning(self):
-        """At the end of the tuning phase, this method
-        removes the first x% of the history so future
-        proposals are not informed by unconverged tuning
-        iterations. Runs only once after the end of tuning,
-        when the self.tuning_end_trigger flag is set to True.
-        """
-        if self.tuning_end_trigger:
-            self.tuning_end_trigger = False
-            return super().stop_tuning()
-        else:
-            return
-
-
-class MLDA(ArrayStepShared):
-    """
-    Multi-Level Delayed Acceptance (MLDA) sampling step that uses coarse
-    approximations of a fine model to construct proposals in multiple levels.
-
-    MLDA creates a hierarchy of MCMC chains. Chains sample from different
-    posteriors that ideally should be approximations of the fine (top-level)
-    posterior and require less computational effort to evaluate their likelihood.
-
-    Each chain runs for a fixed number of iterations (up to subsampling_rate) and
-    then the last sample generated is used as a proposal for the chain in the level
-    above (excluding when variance reduction is used, where a random sample from
-    the generated sequence is used). The bottom-level chain is a MetropolisMLDA
-    or DEMetropolisZMLDA sampler.
-
-    The algorithm achieves higher acceptance rate and effective sample sizes
-    than other samplers if the coarse models are sufficiently good approximations
-    of the fine one.
-
-    Parameters
-    ----------
-    coarse_models : list
-        List of coarse (multi-level) models, where the first model
-        is the coarsest one (level=0) and the last model is the
-        second finest one (level=L-1 where L is the number of levels).
-        Note this list excludes the model passed to the model
-        argument above, which is the finest available.
-    vars : list
-        List of value variables for sampler
-    base_sampler : string
-        Sampler used in the base (coarsest) chain. Can be 'Metropolis' or
-        'DEMetropolisZ'. Defaults to 'DEMetropolisZ'.
-    base_S : standard deviation of base proposal covariance matrix
-        Some measure of variance to parameterize base proposal distribution
-    base_proposal_dist : function
-        Function that returns zero-mean deviates when parameterized with
-        S (and n). Defaults to normal. This is the proposal used in the
-        coarsest (base) chain, i.e. level=0.
-    base_scaling : scalar or array
-        Initial scale factor for base proposal. Defaults to 1 if base_sampler
-        is 'Metropolis' and to 0.001 if base_sampler is 'DEMetropolisZ'.
-    tune : bool
-        Flag for tuning in the base proposal. If base_sampler is 'Metropolis' it
-        should be True or False and defaults to True. Note that
-        this is overidden by the tune parameter in sample(). For example when calling
-        step=MLDA(tune=False, ...) and then sample(step=step, tune=200, ...),
-        tuning will be activated for the first 200 steps. If base_sampler is
-        'DEMetropolisZ', it should be True. For 'DEMetropolisZ', there is a separate
-        argument base_tune_target which allows modifying the type of tuning.
-    base_tune_target: string
-        Defines the type of tuning that is performed when base_sampler is
-        'DEMetropolisZ'. Allowable values are 'lambda, 'scaling' or None and
-        it defaults to 'lambda'.
-    base_tune_interval : int
-        The frequency of tuning for the base proposal. Defaults to 100
-        iterations.
-    base_lamb : float
-        Lambda parameter of the base level DE proposal mechanism. Only applicable when
-        base_sampler is 'DEMetropolisZ'. Defaults to 2.38 / sqrt(2 * ndim)
-    base_tune_drop_fraction : float
-        Fraction of tuning steps that will be removed from the base level samplers
-        history when the tuning ends. Only applicable when base_sampler is
-        'DEMetropolisZ'. Defaults to 0.9 - keeping the last 10% of tuning steps
-        for good mixing while removing 90% of potentially unconverged tuning positions.
-    model : PyMC Model
-        Optional model for sampling step. Defaults to None
-        (taken from context). This model should be the finest of all
-        multilevel models.
-    mode :  string or `Mode` instance.
-        Compilation mode passed to Aesara functions
-    subsampling_rates : integer or list of integers
-        One interger for all levels or a list with one number for each level
-        (excluding the finest level).
-        This is the number of samples generated in level l-1 to
-        propose a sample for level l - applies to all levels excluding the
-        finest level). The length of the list needs to be the same as the
-        length of coarse_models.
-    base_blocked : bool
-        Flag to choose whether base sampler (level=0) is a
-        Compound MetropolisMLDA step (base_blocked=False)
-        or a blocked MetropolisMLDA step (base_blocked=True).
-        Only applicable when base_sampler='Metropolis'.
-    variance_reduction: bool
-        Calculate and store quantities of interest and quantity of interest
-        differences between levels to enable computing a variance-reduced
-        sum of the quantity of interest after sampling. In order to use
-        variance reduction, the user needs to do the following when defining
-        the PyMC model (also demonstrated in the example notebook):
-            - Include a `pm.MutableData()` variable with the name `Q` in the
-            model description of all levels.
-            - Use an Aesara Op to calculate the forward model (or the
-            combination of a forward model and a likelihood). This Op
-            should have a `perform()` method which (in addition to all
-            the other calculations), calculates the quantity of interest
-            and stores it to the variable `Q` of the PyMC model,
-            using the `set_value()` function.
-        When variance_reduction=True, all subchains run for a fixed number
-        of iterations (equal to subsampling_rates) and a random sample is
-        selected from the generated sequence (instead of the last sample
-        which is selected when variance_reduction=False).
-    store_Q_fine: bool
-        Store the values of the quantity of interest from the fine chain.
-    adaptive_error_model : bool
-        When True, MLDA will use the adaptive error model method
-        proposed in [Cui2012]. The method requires the likelihood of
-        the model to be adaptive and a forward model to be defined and
-        fed to the sampler. Thus, it only works when the user does
-        the following (also demonstrated in the example notebook):
-            - Include in the model definition at all levels,
-            the extra variable model_output, which will capture
-            the forward model outputs. Also include in the model
-            definition at all levels except the finest one, the
-            extra variables mu_B and Sigma_B, which will capture
-            the bias between different levels. All these variables
-            should be instantiated using the pm.MutableData method.
-            - Use an Aesara Op to define the forward model (and
-            optionally the likelihood) for all levels. The Op needs
-            to store the result of each forward model calculation
-            to the variable model_output of the PyMC model,
-            using the `set_value()` function.
-            - Define a Multivariate Normal likelihood (either using
-            the standard PyMC API or within an Op) which has mean
-            equal to the forward model output plus mu_B and covariance
-            equal to the model error plus Sigma_B.
-        Given the above, MLDA will capture and iteratively update the
-        bias terms internally for all level pairs and will correct
-        each level so that all levels' forward models aim to estimate
-        the finest level's forward model.
-
-    Examples
-    ----------
-    .. code:: ipython
-        >>> import pymc as pm
-        ... datum = 1
-        ...
-        ... with pm.Model() as coarse_model:
-        ...     x = pm.Normal("x", mu=0, sigma=10)
-        ...     y = pm.Normal("y", mu=x, sigma=1, observed=datum - 0.1)
-        ...
-        ... with pm.Model():
-        ...     x = pm.Normal("x", mu=0, sigma=10)
-        ...     y = pm.Normal("y", mu=x, sigma=1, observed=datum)
-        ...     step_method = pm.MLDA(coarse_models=[coarse_model],
-        ...                           subsampling_rates=5)
-        ...     idata = pm.sample(500, chains=2,
-        ...                       tune=100, step=step_method,
-        ...                       random_seed=123)
-        ...
-        ... az.summary(idata, kind="stats")
-           mean     sd  hdi_3%  hdi_97%
-        x  0.99  0.987  -0.734    2.992
-
-    References
-    ----------
-    .. [Dodwell2019] Dodwell, Tim & Ketelsen, Chris & Scheichl,
-    Robert & Teckentrup, Aretha. (2019).
-    Multilevel Markov Chain Monte Carlo.
-    SIAM Review. 61. 509-545.
-        `link <https://doi.org/10.1137/19M126966X>`__
-    .. [Cui2012] Cui, Tiangang & Fox, Colin & O'Sullivan, Michael.
-    (2012). Adaptive Error Modelling in MCMC Sampling for Large
-    Scale Inverse Problems.
-    """
-
-    name = "mlda"
-
-    # All levels use block sampling,
-    # except level 0 where the user can choose
-    default_blocked = True
-    generates_stats = True
-
-    def __init__(
-        self,
-        coarse_models: List[Model],
-        vars: Optional[list] = None,
-        base_sampler="DEMetropolisZ",
-        base_S: Optional = None,
-        base_proposal_dist: Optional[Type[Proposal]] = None,
-        base_scaling: Optional = None,
-        tune: bool = True,
-        base_tune_target: str = "lambda",
-        base_tune_interval: int = 100,
-        base_lamb: Optional = None,
-        base_tune_drop_fraction: float = 0.9,
-        model: Optional[Model] = None,
-        mode: Optional = None,
-        subsampling_rates: List[int] = 5,
-        base_blocked: bool = False,
-        variance_reduction: bool = False,
-        store_Q_fine: bool = False,
-        adaptive_error_model: bool = False,
-        **kwargs,
-    ) -> None:
-        warnings.warn(
-            "The MLDA stepper will be migrated to the pymc-experimental project!"
-            " See https://github.com/pymc-devs/pymc/issues/5942",
-            DeprecationWarning,
-        )
-
-        # this variable is used to identify MLDA objects which are
-        # not in the finest level (i.e. child MLDA objects)
-        self.is_child = kwargs.get("is_child", False)
-
-        if not isinstance(coarse_models, list):
-            raise ValueError("MLDA step method cannot use coarse_models if it is not a list")
-        if len(coarse_models) == 0:
-            raise ValueError(
-                "MLDA step method was given an empty "
-                "list of coarse models. Give at least "
-                "one coarse model."
-            )
-
-        # assign internal state
-        model = pm.modelcontext(model)
-        initial_values = model.initial_point()
-        self.model = model
-        self.coarse_models = coarse_models
-        self.model_below = self.coarse_models[-1]
-        self.num_levels = len(self.coarse_models) + 1
-
-        # set up variance reduction.
-        self.variance_reduction = variance_reduction
-        self.store_Q_fine = store_Q_fine
-
-        # check that certain requirements hold
-        # for the variance reduction feature to work
-        if self.variance_reduction or self.store_Q_fine:
-            if not hasattr(self.model, "Q"):
-                raise AttributeError(
-                    "Model given to MLDA does not contain"
-                    "variable 'Q'. You need to include"
-                    "the variable in the model definition"
-                    "for variance reduction to work or"
-                    "for storing the fine Q."
-                    "Use pm.MutableData() to define it."
-                )
-            if not isinstance(self.model.Q, TensorSharedVariable):
-                raise TypeError(
-                    "The variable 'Q' in the model definition is not of type "
-                    "'TensorSharedVariable'. Use pm.MutableData() to define the"
-                    "variable."
-                )
-
-        if self.is_child and self.variance_reduction:
-            # this is the subsampling rate applied to the current level
-            # it is stored in the level above and transferred here
-            self.subsampling_rate_above = kwargs.pop("subsampling_rate_above", None)
-
-        # set up adaptive error model
-        self.adaptive_error_model = adaptive_error_model
-
-        # check that certain requirements hold
-        # for the adaptive error model feature to work
-        if self.adaptive_error_model:
-            if not hasattr(self.model_below, "mu_B"):
-                raise AttributeError(
-                    "Model below in hierarchy does not contain"
-                    "variable 'mu_B'. You need to include"
-                    "the variable in the model definition"
-                    "for adaptive error model to work."
-                    "Use pm.MutableData() to define it."
-                )
-            if not hasattr(self.model_below, "Sigma_B"):
-                raise AttributeError(
-                    "Model below in hierarchy does not contain"
-                    "variable 'Sigma_B'. You need to include"
-                    "the variable in the model definition"
-                    "for adaptive error model to work."
-                    "Use pm.MutableData() to define it."
-                )
-            if not (
-                isinstance(self.model_below.mu_B, TensorSharedVariable)
-                and isinstance(self.model_below.Sigma_B, TensorSharedVariable)
-            ):
-                raise TypeError(
-                    "At least one of the variables 'mu_B' and 'Sigma_B' "
-                    "in the definition of the below model is not of type "
-                    "'TensorSharedVariable'. Use pm.MutableData() to define those "
-                    "variables."
-                )
-
-            # this object is used to recursively update the mean and
-            # variance of the bias correction given new differences
-            # between levels
-            self.bias = RecursiveSampleMoments(
-                self.model_below.mu_B.get_value(), self.model_below.Sigma_B.get_value()
-            )
-
-            # this list holds the bias objects from all levels
-            # it is gradually constructed when MLDA objects are
-            # created and then shared between all levels
-            self.bias_all = kwargs.pop("bias_all", None)
-            if self.bias_all is None:
-                self.bias_all = [self.bias]
-            else:
-                self.bias_all.append(self.bias)
-
-            # variables used for adaptive error model
-            self.last_synced_output_diff = None
-            self.adaptation_started = False
-
-        # set up subsampling rates.
-        if isinstance(subsampling_rates, int):
-            self.subsampling_rates = [subsampling_rates] * len(self.coarse_models)
-        else:
-            if len(subsampling_rates) != len(self.coarse_models):
-                raise ValueError(
-                    f"List of subsampling rates needs to have the same "
-                    f"length as list of coarse models but the lengths "
-                    f"were {len(subsampling_rates)}, {len(self.coarse_models)}"
-                )
-            self.subsampling_rates = subsampling_rates
-
-        self.subsampling_rate = self.subsampling_rates[-1]
-        self.subchain_selection = None
-
-        # set up base sampling
-        self.base_sampler = base_sampler
-
-        # VR is not compatible with compound base samplers so an automatic conversion
-        # to a block sampler happens here if
-        if self.variance_reduction and self.base_sampler == "Metropolis" and not base_blocked:
-            warnings.warn(
-                "Variance reduction is not compatible with non-blocked (compound) samplers."
-                "Automatically switching to a blocked Metropolis sampler."
-            )
-            self.base_blocked = True
-        else:
-            self.base_blocked = base_blocked
-
-        self.base_S = base_S
-        self.base_proposal_dist = base_proposal_dist
-
-        if base_scaling is None:
-            if self.base_sampler == "Metropolis":
-                self.base_scaling = 1.0
-            else:
-                self.base_scaling = 0.001
-        else:
-            self.base_scaling = float(base_scaling)
-
-        self.tune = tune
-        if not self.tune and self.base_sampler == "DEMetropolisZ":
-            raise ValueError(
-                f"The argument tune was set to False while using"
-                f" a 'DEMetropolisZ' base sampler. 'DEMetropolisZ' "
-                f" tune needs to be True."
-            )
-
-        self.base_tune_target = base_tune_target
-        self.base_tune_interval = base_tune_interval
-        self.base_lamb = base_lamb
-        self.base_tune_drop_fraction = float(base_tune_drop_fraction)
-        self.base_tuning_stats = None
-
-        self.mode = mode
-
-        # Process model variables
-        if vars is None:
-            vars = model.value_vars
-        else:
-            vars = [model.rvs_to_values.get(var, var) for var in vars]
-        vars = pm.inputvars(vars)
-        self.vars = vars
-        self.var_names = [var.name for var in self.vars]
-
-        self.accepted = 0
-
-        # Construct Aesara function for current-level model likelihood
-        # (for use in acceptance)
-        shared = pm.make_shared_replacements(initial_values, vars, model)
-        self.delta_logp = delta_logp(initial_values, model.logp(), vars, shared)
-
-        # Construct Aesara function for below-level model likelihood
-        # (for use in acceptance)
-        model_below = pm.modelcontext(self.model_below)
-        vars_below = [var for var in model_below.value_vars if var.name in self.var_names]
-        vars_below = pm.inputvars(vars_below)
-        shared_below = pm.make_shared_replacements(initial_values, vars_below, model_below)
-        self.delta_logp_below = delta_logp(
-            initial_values, model_below.logp(), vars_below, shared_below
-        )
-
-        super().__init__(vars, shared)
-
-        # initialise complete step method hierarchy
-        if self.num_levels == 2:
-            with self.model_below:
-                # make sure the correct variables are selected from model_below
-                vars_below = [
-                    var for var in self.model_below.value_vars if var.name in self.var_names
-                ]
-
-                # create kwargs
-                if self.variance_reduction:
-                    base_kwargs = {
-                        "mlda_subsampling_rate_above": self.subsampling_rate,
-                        "mlda_variance_reduction": True,
-                    }
-                else:
-                    base_kwargs = {}
-
-                if self.base_sampler == "Metropolis":
-                    # MetropolisMLDA sampler in base level (level=0), targeting self.model_below
-                    self.step_method_below = pm.MetropolisMLDA(
-                        vars=vars_below,
-                        proposal_dist=self.base_proposal_dist,
-                        S=self.base_S,
-                        scaling=self.base_scaling,
-                        tune=self.tune,
-                        tune_interval=self.base_tune_interval,
-                        model=None,
-                        mode=self.mode,
-                        blocked=self.base_blocked,
-                        **base_kwargs,
-                    )
-                else:
-                    # DEMetropolisZMLDA sampler in base level (level=0), targeting self.model_below
-                    self.step_method_below = pm.DEMetropolisZMLDA(
-                        vars=vars_below,
-                        S=self.base_S,
-                        proposal_dist=self.base_proposal_dist,
-                        lamb=self.base_lamb,
-                        scaling=self.base_scaling,
-                        tune=self.base_tune_target,
-                        tune_interval=self.base_tune_interval,
-                        tune_drop_fraction=self.base_tune_drop_fraction,
-                        model=None,
-                        mode=self.mode,
-                        **base_kwargs,
-                    )
-        else:
-            # drop the last coarse model
-            coarse_models_below = self.coarse_models[:-1]
-            subsampling_rates_below = self.subsampling_rates[:-1]
-
-            with self.model_below:
-                # make sure the correct variables are selected from model_below
-                vars_below = [
-                    var for var in self.model_below.value_vars if var.name in self.var_names
-                ]
-
-                # create kwargs
-                if self.variance_reduction:
-                    mlda_kwargs = {
-                        "is_child": True,
-                        "subsampling_rate_above": self.subsampling_rate,
-                    }
-                else:
-                    mlda_kwargs = {"is_child": True}
-                if self.adaptive_error_model:
-                    mlda_kwargs = {**mlda_kwargs, **{"bias_all": self.bias_all}}
-
-                # MLDA sampler in some intermediate level, targeting self.model_below
-                self.step_method_below = pm.MLDA(
-                    vars=vars_below,
-                    base_S=self.base_S,
-                    base_sampler=self.base_sampler,
-                    base_proposal_dist=self.base_proposal_dist,
-                    base_scaling=self.base_scaling,
-                    tune=self.tune,
-                    base_tune_target=self.base_tune_target,
-                    base_tune_interval=self.base_tune_interval,
-                    base_lamb=self.base_lamb,
-                    base_tune_drop_fraction=self.base_tune_drop_fraction,
-                    model=None,
-                    mode=self.mode,
-                    subsampling_rates=subsampling_rates_below,
-                    coarse_models=coarse_models_below,
-                    base_blocked=self.base_blocked,
-                    variance_reduction=self.variance_reduction,
-                    store_Q_fine=False,
-                    adaptive_error_model=self.adaptive_error_model,
-                    **mlda_kwargs,
-                )
-
-        # instantiate the recursive DA proposal.
-        # this is the main proposal used for
-        # all levels (Recursive Delayed Acceptance)
-        # (except for level 0 where the step method is MetropolisMLDA
-        # or DEMetropolisZMLDA - not MLDA)
-        self.proposal_dist = RecursiveDAProposal(
-            self.step_method_below, self.model_below, self.tune, self.subsampling_rate
-        )
-
-        # set up data types of stats.
-        if isinstance(self.step_method_below, MLDA):
-            # get the stat types from the level below if that level is MLDA
-            self.stats_dtypes = self.step_method_below.stats_dtypes
-
-        else:
-            # otherwise, set it up from scratch.
-            self.stats_dtypes = [{"accept": np.float64, "accepted": bool, "tune": bool}]
-
-            if isinstance(self.step_method_below, MetropolisMLDA):
-                self.stats_dtypes.append({"base_scaling": np.float64})
-            elif isinstance(self.step_method_below, DEMetropolisZMLDA):
-                self.stats_dtypes.append({"base_scaling": np.float64, "base_lambda": np.float64})
-            elif isinstance(self.step_method_below, CompoundStep):
-                for method in self.step_method_below.methods:
-                    if isinstance(method, MetropolisMLDA):
-                        self.stats_dtypes.append({"base_scaling": np.float64})
-                    elif isinstance(method, DEMetropolisZMLDA):
-                        self.stats_dtypes.append(
-                            {"base_scaling": np.float64, "base_lambda": np.float64}
-                        )
-
-        # initialise necessary variables for doing variance reduction
-        if self.variance_reduction:
-            self.sub_counter = 0
-            self.Q_diff = []
-            if self.is_child:
-                self.Q_reg = [np.nan] * self.subsampling_rate_above
-            if self.num_levels == 2:
-                self.Q_base_full = []
-            if not self.is_child:
-                for level in range(self.num_levels - 1, 0, -1):
-                    self.stats_dtypes[0][f"Q_{level}_{level - 1}"] = object
-                self.stats_dtypes[0]["Q_0"] = object
-
-        # initialise necessary variables for doing variance reduction or storing fine Q
-        if self.variance_reduction or self.store_Q_fine:
-            self.Q_last = np.nan
-            self.Q_diff_last = np.nan
-        if self.store_Q_fine and not self.is_child:
-            self.stats_dtypes[0][f"Q_{self.num_levels - 1}"] = object
-
-    def astep(self, q0: RaveledVars) -> Tuple[RaveledVars, List[Dict[str, Any]]]:
-        """One MLDA step, given current sample q0"""
-        # Check if the tuning flag has been changed and if yes,
-        # change the proposal's tuning flag and reset self.accepted
-        # This is triggered by _iter_sample while the highest-level MLDA step
-        # method is running. It then propagates to all levels.
-        if self.proposal_dist.tune != self.tune:
-            self.proposal_dist.tune = self.tune
-            # set tune in sub-methods of compound stepper explicitly because
-            # it is not set within sample.py (only the CompoundStep's tune flag is)
-            if isinstance(self.step_method_below, CompoundStep):
-                for method in self.step_method_below.methods:
-                    method.tune = self.tune
-            self.accepted = 0
-
-        # Set subchain_selection (which sample from the coarse chain
-        # is passed as a proposal to the fine chain). If variance
-        # reduction is used, a random sample is selected as proposal.
-        # If variance reduction is not used, the last sample is
-        # selected as proposal.
-        if self.variance_reduction:
-            self.subchain_selection = np.random.randint(0, self.subsampling_rate)
-        else:
-            self.subchain_selection = self.subsampling_rate - 1
-        self.proposal_dist.subchain_selection = self.subchain_selection
-
-        # Call the recursive DA proposal to get proposed sample
-        # and convert dict -> numpy array
-        q = self.proposal_dist(q0)
-
-        # Evaluate MLDA acceptance log-ratio
-        # If proposed sample from lower levels is the same as current one,
-        # do not calculate likelihood, just set accept to 0.0
-        if (q.data == q0.data).all():
-            accept = np.float64(0.0)
-            skipped_logp = True
-        else:
-            # NB! The order and sign of the first term are swapped compared
-            # to the convention to make sure the proposal is evaluated last.
-            accept = -self.delta_logp(q0.data, q.data) + self.delta_logp_below(q0.data, q.data)
-            skipped_logp = False
-
-        # Accept/reject sample - next sample is stored in q_new
-        q_new, accepted = metrop_select(accept, q, q0)
-        if skipped_logp:
-            accepted = False
-
-        # if sample is accepted, update self.Q_last with the sample's Q value
-        # runs only for VR or when store_Q_fine is True
-        if self.variance_reduction or self.store_Q_fine:
-            if accepted and not skipped_logp:
-                self.Q_last = self.model.Q.get_value()
-
-        # Variance reduction
-        if self.variance_reduction:
-            self.update_vr_variables(accepted, skipped_logp)
-
-        # Adaptive error model - runs only during tuning.
-        if self.tune and self.adaptive_error_model:
-            self.update_error_estimate(accepted, skipped_logp)
-
-        # Update acceptance counter
-        self.accepted += accepted
-
-        stats = {"tune": self.tune, "accept": np.exp(accept), "accepted": accepted}
-
-        # Save the VR statistics to the stats dictionary (only happens in the
-        # top MLDA level)
-        if (self.variance_reduction or self.store_Q_fine) and not self.is_child:
-            q_stats = {}
-            if self.variance_reduction:
-                m = self
-                for level in range(self.num_levels - 1, 0, -1):
-                    # save the Q differences for this level and iteration
-                    q_stats[f"Q_{level}_{level - 1}"] = np.array(m.Q_diff)
-                    # this makes sure Q_diff is reset for
-                    # the next iteration
-                    m.Q_diff = []
-                    if level == 1:
-                        break
-                    m = m.step_method_below
-                q_stats["Q_0"] = np.array(m.Q_base_full)
-                m.Q_base_full = []
-            if self.store_Q_fine:
-                q_stats["Q_" + str(self.num_levels - 1)] = np.array(self.Q_last)
-            stats = {**stats, **q_stats}
-
-        # Capture the base tuning stats from the level below.
-        self.base_tuning_stats = []
-
-        if isinstance(self.step_method_below, MLDA):
-            self.base_tuning_stats = self.step_method_below.base_tuning_stats
-        elif isinstance(self.step_method_below, MetropolisMLDA):
-            self.base_tuning_stats.append({"base_scaling": np.mean(self.step_method_below.scaling)})
-        elif isinstance(self.step_method_below, DEMetropolisZMLDA):
-            self.base_tuning_stats.append(
-                {
-                    "base_scaling": np.mean(self.step_method_below.scaling),
-                    "base_lambda": self.step_method_below.lamb,
-                }
-            )
-        elif isinstance(self.step_method_below, CompoundStep):
-            # Below method is CompoundStep
-            for method in self.step_method_below.methods:
-                if isinstance(method, MetropolisMLDA):
-                    self.base_tuning_stats.append({"base_scaling": np.mean(method.scaling)})
-                elif isinstance(method, DEMetropolisZMLDA):
-                    self.base_tuning_stats.append(
-                        {"base_scaling": np.mean(method.scaling), "base_lambda": method.lamb}
-                    )
-
-        return q_new, [stats] + self.base_tuning_stats
-
-    def update_vr_variables(self, accepted, skipped_logp):
-        """Updates all the variables necessary for VR to work.
-
-        Each level has a Q_last and Q_diff_last register which store
-        the Q of the last accepted MCMC sample and the difference
-        between the Q of the last accepted sample in this level and
-        the Q of the last sample in the level below.
-
-        These registers are updated here so that they can be exported later."""
-
-        # if this MLDA is not at the finest level, store Q_last in a
-        # register Q_reg and increase sub_counter (until you reach
-        # the subsampling rate, at which point you make it zero).
-        # Q_reg will later be used by the level above to calculate differences
-        if self.is_child:
-            if self.sub_counter == self.subsampling_rate_above:
-                self.sub_counter = 0
-            self.Q_reg[self.sub_counter] = self.Q_last
-            self.sub_counter += 1
-
-        # if MLDA is in the level above the base level, extract the
-        # latest set of Q values from Q_reg in the base level
-        # and add them to Q_base_full (which stores all the history of
-        # Q values from the base level)
-        if self.num_levels == 2:
-            self.Q_base_full.extend(self.step_method_below.Q_reg)
-
-        # if the sample is accepted, update Q_diff_last with the latest
-        # difference between the last Q of this level and the Q of the
-        # proposed (selected) sample from the level below.
-        # If sample is not accepted, just keep the latest accepted Q_diff
-        if accepted and not skipped_logp:
-            self.Q_diff_last = self.Q_last - self.step_method_below.Q_reg[self.subchain_selection]
-        # Add the last accepted Q_diff to the list
-        self.Q_diff.append(self.Q_diff_last)
-
-    def update_error_estimate(self, accepted, skipped_logp):
-        """Updates the adaptive error model estimate with
-        the latest accepted forward model output difference. Also
-        updates the model variables mu_B and Sigma_B.
-
-        The current level estimates and stores the error
-        model between the current level and the level below."""
-
-        # only save errors when a sample is accepted (excluding skipped_logp)
-        if accepted and not skipped_logp:
-            # this is the error (i.e. forward model output difference)
-            # between the current level's model and the model in the level below
-            self.last_synced_output_diff = (
-                self.model.model_output.get_value() - self.model_below.model_output.get_value()
-            )
-            self.adaptation_started = True
-
-        if self.adaptation_started:
-            # update the internal recursive bias estimator with the last saved error
-            self.bias.update(self.last_synced_output_diff)
-            # Update the model variables in the level below the current one.
-            # Each level has its own bias correction (i.e. bias object) that
-            # estimates the error between that level and the one below.
-            # The model variables mu_B and Signa_B of a level are the
-            # sum of the bias corrections of all levels below and including
-            # that level. This sum is updated here.
-            with self.model_below:
-                pm.set_data(
-                    {
-                        "mu_B": sum(
-                            bias.get_mu()
-                            for bias in self.bias_all[: len(self.bias_all) - self.num_levels + 2]
-                        )
-                    }
-                )
-                pm.set_data(
-                    {
-                        "Sigma_B": sum(
-                            bias.get_sigma()
-                            for bias in self.bias_all[: len(self.bias_all) - self.num_levels + 2]
-                        )
-                    }
-                )
-
-    @staticmethod
-    def competence(var, has_grad):
-        """Return MLDA competence for given var/has_grad. MLDA currently works
-        only with continuous variables."""
-        if var.dtype in pm.discrete_types:
-            return Competence.INCOMPATIBLE
-        return Competence.COMPATIBLE
-
-
-class RecursiveSampleMoments:
-    """
-    Iteratively constructs a sample mean
-    and covariance, given input samples.
-
-    Used to capture an estimate of the mean
-    and covariance of the bias of an MLDA
-    coarse model.
-    """
-
-    def __init__(self, mu_0, sigma_0, t=1):
-        self.mu = mu_0
-        self.sigma = sigma_0
-        self.t = t
-
-    def __call__(self):
-        return self.mu, self.sigma
-
-    def get_mu(self):
-        """Returns the current mu value"""
-        return self.mu
-
-    def get_sigma(self):
-        """Returns the current covariance value"""
-        return self.sigma
-
-    def update(self, x):
-        """Updates the mean and covariance given a
-        new sample x"""
-        mu_previous = self.mu.copy()
-
-        self.mu = (1 / (self.t + 1)) * (self.t * mu_previous + x)
-
-        self.sigma = (self.t - 1) / self.t * self.sigma + 1 / self.t * (
-            self.t * np.outer(mu_previous, mu_previous)
-            - (self.t + 1) * np.outer(self.mu, self.mu)
-            + np.outer(x, x)
-        )
-
-        self.t += 1
-
-
-def extract_Q_estimate(trace, levels):
-    """
-    Returns expectation and standard error of quantity of interest,
-    given a trace and the number of levels in the multilevel model.
-    It makes use of the collapsing sum formula. Only applicable when
-    MLDA with variance reduction has been used for sampling.
-    """
-
-    Q_0_raw = trace.get_sampler_stats("Q_0").squeeze()
-    Q_0 = np.concatenate(Q_0_raw)[None, ::]
-    ess_Q_0 = az.ess(np.array(Q_0, np.float64))
-    Q_0_var = Q_0.var() / ess_Q_0
-
-    Q_diff_means = []
-    Q_diff_vars = []
-    for l in range(1, levels):
-        Q_diff_raw = trace.get_sampler_stats(f"Q_{l}_{l-1}").squeeze()
-        Q_diff = np.hstack(Q_diff_raw)[None, ::]
-        ess_diff = az.ess(np.array(Q_diff, np.float64))
-        Q_diff_means.append(Q_diff.mean())
-        Q_diff_vars.append(Q_diff.var() / ess_diff)
-
-    Q_mean = Q_0.mean() + sum(Q_diff_means)
-    Q_se = np.sqrt(Q_0_var + sum(Q_diff_vars))
-    return Q_mean, Q_se
-
-
-def subsample(
-    draws=1,
-    step=None,
-    start=None,
-    trace=None,
-    tune=0,
-    model=None,
-):
-    """
-    A stripped down version of sample(), which is called only
-    by the RecursiveDAProposal (which is the proposal used in the MLDA
-    stepper). RecursiveDAProposal only requires a small set of the input
-    parameters and checks normally performed by sample(), and this
-    function thus skips some of the code in sampler(). It directly calls
-    _iter_sample(), rather than sample_many(). The result is a reduced
-    overhead when running multiple levels in MLDA.
-    """
-
-    model = pm.modelcontext(model)
-    chain = 0
-    random_seed = np.random.randint(2**30)
-    callback = None
-
-    draws += tune
-
-    sampling = pm.sampling._iter_sample(
-        draws, step, start, trace, chain, tune, model, random_seed, callback
-    )
-
-    try:
-        for it, (trace, _) in enumerate(sampling):
-            pass
-    except KeyboardInterrupt:
-        pass
-
-    return trace
-
-
-# Available proposal distributions for MLDA
-
-
-class RecursiveDAProposal(Proposal):
-    """
-    Recursive Delayed Acceptance proposal to be used with MLDA step sampler.
-    Recursively calls an MLDA sampler if level > 0 and calls MetropolisMLDA or
-    DEMetropolisZMLDA sampler if level = 0. The sampler generates
-    self.subsampling_rate samples and returns the sample with index
-    self.subchain_selection to be used as a proposal.
-    Results in a hierarchy of chains each of which is used to propose
-    samples to the chain above.
-    """
-
-    def __init__(
-        self,
-        step_method_below: Union[MLDA, MetropolisMLDA, DEMetropolisZMLDA, CompoundStep],
-        model_below: Model,
-        tune: bool,
-        subsampling_rate: int,
-    ) -> None:
-
-        self.step_method_below = step_method_below
-        self.model_below = model_below
-        self.tune = tune
-        self.subsampling_rate = subsampling_rate
-        self.subchain_selection = None
-        self.tuning_end_trigger = True
-
-    def __call__(self, q0: RaveledVars) -> RaveledVars:
-        """Returns proposed sample given the current sample
-        in dictionary form (q0_dict)."""
-
-        # Logging is reduced to avoid extensive console output
-        # during multiple recursive calls of subsample()
-        _log = logging.getLogger("pymc")
-        _log.setLevel(logging.ERROR)
-
-        # Convert current sample from RaveledVars ->
-        # dict before feeding to subsample.
-        q0_dict = DictToArrayBijection.rmap(q0)
-
-        with self.model_below:
-            # Check if the tuning flag has been set to False
-            # in which case tuning is stopped. The flag is set
-            # to False (by MLDA's astep) when the burn-in
-            # iterations of the highest-level MLDA sampler run out.
-            # The change propagates to all levels.
-
-            if self.tune:
-                # Subsample in tuning mode
-                trace = subsample(
-                    draws=0,
-                    step=self.step_method_below,
-                    start=q0_dict,
-                    tune=self.subsampling_rate,
-                )
-            else:
-                # Subsample in normal mode without tuning
-                # If DEMetropolisZMLDA is the base sampler a flag is raised to
-                # make sure that history is edited after tuning ends
-                if self.tuning_end_trigger:
-                    if isinstance(self.step_method_below, DEMetropolisZMLDA):
-                        self.step_method_below.tuning_end_trigger = True
-                    self.tuning_end_trigger = False
-
-                trace = subsample(
-                    draws=self.subsampling_rate,
-                    step=self.step_method_below,
-                    start=q0_dict,
-                    tune=0,
-                )
-
-        # set logging back to normal
-        _log.setLevel(logging.NOTSET)
-
-        # return sample with index self.subchain_selection from the generated
-        # sequence of length self.subsampling_rate. The index is set within
-        # MLDA's astep() function
-        q_dict = trace.point(self.subchain_selection)
-
-        # Make sure output dict is ordered the same way as the input dict.
-        q_dict = Point(
-            {key: q_dict[key] for key in q0_dict.keys()},
-            model=self.model_below,
-            filter_model_vars=True,
-        )
-
-        return DictToArrayBijection.map(q_dict)
diff --git a/pymc/util.py b/pymc/util.py
index 657b9a8f5..8fc6f3cf7 100644
--- a/pymc/util.py
+++ b/pymc/util.py
@@ -14,7 +14,7 @@
 
 import functools
 
-from typing import Dict, List, Tuple, Union, cast
+from typing import Dict, Hashable, List, Tuple, Union, cast
 
 import arviz
 import cloudpickle
@@ -232,15 +232,18 @@ def biwrap(wrapper):
 
 
 def dataset_to_point_list(ds: xarray.Dataset) -> List[Dict[str, np.ndarray]]:
+    # All keys of the dataset must be a str
+    for vn in ds.keys():
+        if not isinstance(vn, str):
+            raise ValueError(f"Variable names must be str, but dataset key {vn} is a {type(vn)}.")
     # make dicts
-    points: List[Dict[str, np.ndarray]] = []
-    vn: str
+    points: List[Dict[Hashable, np.ndarray]] = []
     da: "xarray.DataArray"
     for c in ds.chain:
         for d in ds.draw:
             points.append({vn: da.sel(chain=c, draw=d).values for vn, da in ds.items()})
     # use the list of points
-    return points
+    return cast(List[Dict[str, np.ndarray]], points)
 
 
 def chains_and_samples(data: Union[xarray.Dataset, arviz.InferenceData]) -> Tuple[int, int]:
diff --git a/requirements-dev.txt b/requirements-dev.txt
index b48830961..d39708643 100644
--- a/requirements-dev.txt
+++ b/requirements-dev.txt
@@ -2,7 +2,7 @@
 # See that file for comments about the need/usage of each dependency.
 
 aeppl==0.0.32
-aesara==2.7.7
+aesara==2.7.8
 arviz>=0.12.0
 cachetools>=4.2.1
 cloudpickle
diff --git a/requirements.txt b/requirements.txt
index 479ca4449..73a70fb98 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -1,5 +1,5 @@
 aeppl==0.0.32
-aesara==2.7.7
+aesara==2.7.8
 arviz>=0.12.0
 cachetools>=4.2.1
 cloudpickle
diff --git a/scripts/run_mypy.py b/scripts/run_mypy.py
index 55c5ef764..8539b3ff1 100644
--- a/scripts/run_mypy.py
+++ b/scripts/run_mypy.py
@@ -101,7 +101,7 @@ def mypy_to_pandas(input_lines: Iterator[str]) -> pandas.DataFrame:
         "file": [],
         "line": [],
         "type": [],
-        "section": [],
+        "errorcode": [],
         "message": [],
     }
     for line in input_lines:
@@ -120,7 +120,7 @@ def mypy_to_pandas(input_lines: Iterator[str]) -> pandas.DataFrame:
             data["file"].append(file)
             data["line"].append(lineno)
             data["type"].append(message_type)
-            data["section"].append(current_section)
+            data["errorcode"].append(current_section)
             data["message"].append(message)
         except Exception as ex:
             print(elems)
@@ -154,21 +154,29 @@ def check_no_unexpected_results(mypy_lines: Iterator[str]):
     if not unexpected_failing:
         print(f"{len(passing)}/{len(all_files)} files pass as expected.")
     else:
-        print(f"{len(unexpected_failing)} files unexpectedly failed:")
+        print("!!!!!!!!!")
+        print(f"{len(unexpected_failing)} files unexpectedly failed.")
         print("\n".join(sorted(map(str, unexpected_failing))))
+        print(
+            "These files did not fail before, so please check the above output"
+            f" for errors in {unexpected_failing} and fix them."
+        )
+        print("!!!!!!!!!")
         sys.exit(1)
 
     if unexpected_passing == {"pymc/sampling_jax.py"}:
         print("Letting you know that 'pymc/sampling_jax.py' unexpectedly passed.")
         print("But this file is known to sometimes pass and sometimes not.")
-        print("Unless tried to resolve problems in sampling_jax.py just ignore this message.")
+        print("Unless you tried to resolve problems in sampling_jax.py just ignore this message.")
     elif unexpected_passing:
+        print("!!!!!!!!!")
         print(f"{len(unexpected_passing)} files unexpectedly passed the type checks:")
         print("\n".join(sorted(map(str, unexpected_passing))))
         print("This is good news! Go to scripts/run_mypy.py and add them to the list.")
         if all_files.issubset(passing):
             print("WOW! All files are passing the mypy type checks!")
             print("scripts\\run_mypy.py may no longer be needed.")
+        print("!!!!!!!!!")
         sys.exit(1)
     return
 
@@ -176,21 +184,37 @@ def check_no_unexpected_results(mypy_lines: Iterator[str]):
 if __name__ == "__main__":
     # Enforce PEP 561 for some important dependencies that
     # have relevant type hints but don't tell that to mypy.
-    enforce_pep561("aesara")
     enforce_pep561("aeppl")
 
-    parser = argparse.ArgumentParser(description="Process some integers.")
+    parser = argparse.ArgumentParser(description="Run mypy type checks on PyMC codebase.")
     parser.add_argument(
         "--verbose", action="count", default=0, help="Pass this to print mypy output."
     )
+    parser.add_argument(
+        "--groupby",
+        default="file",
+        help="How to group verbose output. One of {file|errorcode|message}.",
+    )
     args, _ = parser.parse_known_args()
 
     cp = subprocess.run(
-        ["mypy", "--exclude", "pymc/tests", "pymc"],
+        ["mypy", "--show-error-codes", "--exclude", "pymc/tests", "pymc"],
         capture_output=True,
     )
     output = cp.stdout.decode()
     if args.verbose:
-        print(output)
+        df = mypy_to_pandas(output.split("\n"))
+        for section, sdf in df.reset_index().groupby(args.groupby):
+            print(f"\n\n[{section}]")
+            for row in sdf.itertuples():
+                print(f"{row.file}:{row.line}: {row.type}: {row.message}")
+        print()
+    else:
+        print(
+            "Mypy output hidden."
+            " Run `python run_mypy.py --verbose` to see the full output,"
+            " or `python run_mypy.py --help` for other options."
+        )
+
     check_no_unexpected_results(output.split("\n"))
     sys.exit(0)
diff --git a/setup.py b/setup.py
index 755f5dc5b..25870bfd6 100755
--- a/setup.py
+++ b/setup.py
@@ -29,6 +29,7 @@ classifiers = [
     "Development Status :: 5 - Production/Stable",
     "Programming Language :: Python",
     "Programming Language :: Python :: 3",
+    "Programming Language :: Python :: 3.7",
     "Programming Language :: Python :: 3.8",
     "Programming Language :: Python :: 3.9",
     "Programming Language :: Python :: 3.10",
@@ -84,7 +85,7 @@ if __name__ == "__main__":
         # package_data={'docs': ['*']},
         include_package_data=True,
         classifiers=classifiers,
-        python_requires=">=3.8",
+        python_requires=">=3.7",
         install_requires=install_reqs,
         tests_require=test_reqs,
     )


Test patch:
diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml
index f479dcc2d..2958e1010 100644
--- a/.github/workflows/tests.yml
+++ b/.github/workflows/tests.yml
@@ -417,3 +417,108 @@ jobs:
           env_vars: TEST_SUBSET
           name: ${{ matrix.os }} ${{ matrix.floatx }}
           fail_ci_if_error: false
+  ubuntu_37:
+    strategy:
+      matrix:
+        os: [ubuntu-20.04]
+        floatx: [float64]
+        python-version: ["3.7"]
+        test-subset:
+          - |
+            pymc/tests/test_util.py
+            pymc/tests/test_logprob.py
+            pymc/tests/test_aesaraf.py
+            pymc/tests/test_math.py
+            pymc/tests/test_posdef_sym.py
+            pymc/tests/test_ndarray_backend.py
+            pymc/tests/test_hmc.py
+            pymc/tests/test_func_utils.py
+            pymc/tests/test_shape_handling.py
+            pymc/tests/test_starting.py
+            pymc/tests/test_mixture.py
+
+          - |
+            pymc/tests/test_distributions.py
+
+          - |
+            pymc/tests/test_tuning.py
+            pymc/tests/test_shared.py
+            pymc/tests/test_types.py
+            pymc/tests/test_distributions_moments.py
+
+          - |
+            pymc/tests/test_modelcontext.py
+            pymc/tests/test_dist_math.py
+            pymc/tests/test_minibatches.py
+            pymc/tests/test_pickling.py
+            pymc/tests/test_transforms.py
+            pymc/tests/test_parallel_sampling.py
+
+          - |
+            pymc/tests/test_distributions_random.py
+            pymc/tests/test_distributions_timeseries.py
+            pymc/tests/test_gp.py
+            pymc/tests/test_model.py
+            pymc/tests/test_model_graph.py
+            pymc/tests/test_ode.py
+            pymc/tests/test_profile.py
+            pymc/tests/test_quadpotential.py
+
+      fail-fast: false
+    runs-on: ${{ matrix.os }}
+    env:
+      TEST_SUBSET: ${{ matrix.test-subset }}
+      AESARA_FLAGS: floatX=${{ matrix.floatx }},gcc__cxxflags='-march=native'
+    defaults:
+      run:
+        shell: bash -l {0}
+    steps:
+      - uses: actions/checkout@v2
+      - name: Cache conda
+        uses: actions/cache@v3
+        env:
+          # Increase this value to reset cache if environment-test.yml has not changed
+          CACHE_NUMBER: 0
+        with:
+          path: ~/conda_pkgs_dir
+          key: ${{ runner.os }}-py${{matrix.python-version}}-conda-${{ env.CACHE_NUMBER }}-${{
+            hashFiles('conda-envs/environment-test.yml') }}
+      - name: Cache multiple paths
+        uses: actions/cache@v3
+        env:
+          # Increase this value to reset cache if requirements.txt has not changed
+          CACHE_NUMBER: 0
+        with:
+          path: |
+            ~/.cache/pip
+            $RUNNER_TOOL_CACHE/Python/*
+            ~\AppData\Local\pip\Cache
+          key: ${{ runner.os }}-build-${{ matrix.python-version }}-${{ env.CACHE_NUMBER }}-${{
+            hashFiles('requirements.txt') }}
+      - uses: conda-incubator/setup-miniconda@v2
+        with:
+          miniforge-variant: Mambaforge
+          miniforge-version: latest
+          mamba-version: "*"
+          activate-environment: pymc-test
+          channel-priority: strict
+          environment-file: conda-envs/environment-test.yml
+          python-version: ${{matrix.python-version}}
+          use-mamba: true
+          use-only-tar-bz2: true # IMPORTANT: This needs to be set for caching to work properly!
+      - name: Install-pymc
+        run: |
+          conda activate pymc-test
+          pip install -e .
+          pip install --pre -U polyagamma
+          python --version
+      - name: Run tests
+        run: |
+          conda activate pymc-test
+          python -m pytest -vv --cov=pymc --cov-report=xml --cov-report term --durations=50 $TEST_SUBSET
+      - name: Upload coverage to Codecov
+        uses: codecov/codecov-action@v2
+        with:
+          env_vars: TEST_SUBSET
+          name: ${{ matrix.os }} ${{ matrix.floatx }}
+          fail_ci_if_error: false
diff --git a/conda-envs/environment-test.yml b/conda-envs/environment-test.yml
index ccac4652f..efb3477f7 100644
--- a/conda-envs/environment-test.yml
+++ b/conda-envs/environment-test.yml
@@ -6,7 +6,7 @@ channels:
 dependencies:
 # Base dependencies
 - aeppl=0.0.32
-- aesara=2.7.7
+- aesara=2.7.8
 - arviz>=0.12.0
 - blas
 - cachetools>=4.2.1
diff --git a/conda-envs/windows-environment-test.yml b/conda-envs/windows-environment-test.yml
index a373ec607..164a29f94 100644
--- a/conda-envs/windows-environment-test.yml
+++ b/conda-envs/windows-environment-test.yml
@@ -6,7 +6,7 @@ channels:
 dependencies:
 # Base dependencies (see install guide for Windows)
 - aeppl=0.0.32
-- aesara=2.7.7
+- aesara=2.7.8
 - arviz>=0.12.0
 - blas
 - cachetools>=4.2.1
diff --git a/pymc/tests/test_distributions_moments.py b/pymc/tests/test_distributions_moments.py
index 45d9db503..75c29849b 100644
--- a/pymc/tests/test_distributions_moments.py
+++ b/pymc/tests/test_distributions_moments.py
@@ -892,19 +892,32 @@ def test_categorical_moment(p, size, expected):
             np.array([-4, -1, 3, 9, 19]),
             np.array([0.1, 0.15, 0.2, 0.25, 0.3]),
             None,
-            1.5458937198067635,
+            9.34782609,
         ),
         (
             np.array([-22, -4, 0, 8, 13]),
             np.tile(1 / 5, 5),
             (5, 3),
-            np.full((5, 3), -0.14285714285714296),
+            np.full((5, 3), -4.5),
         ),
         (
             np.arange(-100, 10),
             np.arange(1, 111) / 6105,
             (2, 5, 3),
-            np.full((2, 5, 3), -27.584097859327223),
+            np.full((2, 5, 3), -27.65765766),
+        ),
+        (
+            # from https://github.com/pymc-devs/pymc/issues/5959
+            np.linspace(0, 10, 10),
+            st.norm.pdf(np.linspace(0, 10, 10), loc=2.5, scale=1),
+            None,
+            2.5270134,
+        ),
+        (
+            np.linspace(0, 10, 100),
+            st.norm.pdf(np.linspace(0, 10, 100), loc=2.5, scale=1),
+            None,
+            2.51771721,
         ),
     ],
 )
diff --git a/pymc/tests/test_step.py b/pymc/tests/test_step.py
index 49abd01e8..8cf809069 100644
--- a/pymc/tests/test_step.py
+++ b/pymc/tests/test_step.py
@@ -24,12 +24,10 @@ import numpy.testing as npt
 import pytest
 
 from aesara.compile.ops import as_op
-from aesara.graph.op import Op
 
 import pymc as pm
 
 from pymc.aesaraf import floatX
-from pymc.data import Data, MutableData
 from pymc.distributions import (
     Bernoulli,
     Beta,
@@ -38,14 +36,12 @@ from pymc.distributions import (
     Dirichlet,
     HalfNormal,
     Multinomial,
-    MvNormal,
     Normal,
 )
 from pymc.exceptions import SamplingError
-from pymc.model import Model, Potential, set_data
+from pymc.model import Model
 from pymc.sampling import assign_step_methods, sample
 from pymc.step_methods import (
-    MLDA,
     NUTS,
     BinaryGibbsMetropolis,
     CategoricalGibbsMetropolis,
@@ -56,18 +52,14 @@ from pymc.step_methods import (
     Metropolis,
     MultivariateNormalProposal,
     NormalProposal,
-    RecursiveDAProposal,
     Slice,
-    UniformProposal,
 )
-from pymc.step_methods.mlda import extract_Q_estimate
 from pymc.tests.checks import close_to
 from pymc.tests.helpers import fast_unstable_sampling_mode
 from pymc.tests.models import (
     mv_simple,
     mv_simple_coarse,
     mv_simple_discrete,
-    mv_simple_very_coarse,
     simple_2model_continuous,
     simple_categorical,
 )
@@ -111,15 +103,6 @@ class TestStepMethods:
                 ),
                 1000,
             ),
-            # MLDA takes 1/2 of the total test time!
-            (
-                lambda C, model_coarse: MLDA(
-                    coarse_models=[model_coarse],
-                    base_S=C,
-                    base_proposal_dist=MultivariateNormalProposal,
-                ),
-                1000,
-            ),
             (lambda *_: Slice(), 2000),
             (lambda *_: Slice(blocked=True), 2000),
             (
@@ -662,706 +645,6 @@ class TestNutsCheckTrace:
         assert (trace.model_logp == model_logp_).all()
 
 
-class TestMLDA:
-    steppers = [MLDA]
-
-    def test_proposal_and_base_proposal_choice(self):
-        """Test that proposal_dist and base_proposal_dist are set as
-        expected by MLDA"""
-        _, model, _ = mv_simple()
-        _, model_coarse, _ = mv_simple_coarse()
-        with model:
-            sampler = MLDA(coarse_models=[model_coarse], base_sampler="Metropolis")
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, NormalProposal)
-
-            sampler = MLDA(coarse_models=[model_coarse])
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, UniformProposal)
-
-            initial_point = model.initial_point()
-            initial_point_size = sum(initial_point[n.name].size for n in model.value_vars)
-            s = np.ones(initial_point_size)
-            sampler = MLDA(coarse_models=[model_coarse], base_sampler="Metropolis", base_S=s)
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, NormalProposal)
-
-            sampler = MLDA(coarse_models=[model_coarse], base_S=s)
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, UniformProposal)
-
-            s = np.diag(s)
-            sampler = MLDA(coarse_models=[model_coarse], base_sampler="Metropolis", base_S=s)
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, MultivariateNormalProposal)
-
-            sampler = MLDA(coarse_models=[model_coarse], base_S=s)
-            assert isinstance(sampler.proposal_dist, RecursiveDAProposal)
-            assert sampler.base_proposal_dist is None
-            assert isinstance(sampler.step_method_below.proposal_dist, UniformProposal)
-
-            s[0, 0] = -s[0, 0]
-            with pytest.raises(np.linalg.LinAlgError):
-                MLDA(coarse_models=[model_coarse], base_sampler="Metropolis", base_S=s)
-
-    def test_step_methods_in_each_level(self):
-        """Test that MLDA creates the correct hierarchy of step methods when no
-        coarse models are passed and when two coarse models are passed."""
-        _, model, _ = mv_simple()
-        _, model_coarse, _ = mv_simple_coarse()
-        _, model_very_coarse, _ = mv_simple_very_coarse()
-        with model:
-            initial_point = model.initial_point()
-            initial_point_size = sum(initial_point[n.name].size for n in model.value_vars)
-            s = np.ones(initial_point_size) + 2.0
-            sampler = MLDA(
-                coarse_models=[model_very_coarse, model_coarse],
-                base_S=s,
-                base_sampler="Metropolis",
-            )
-            assert isinstance(sampler.step_method_below, MLDA)
-            assert isinstance(sampler.step_method_below.step_method_below, Metropolis)
-            assert np.all(sampler.step_method_below.step_method_below.proposal_dist.s == s)
-
-            sampler = MLDA(coarse_models=[model_very_coarse, model_coarse], base_S=s)
-            assert isinstance(sampler.step_method_below, MLDA)
-            assert isinstance(sampler.step_method_below.step_method_below, DEMetropolisZ)
-            assert np.all(sampler.step_method_below.step_method_below.proposal_dist.s == s)
-
-    def test_exceptions_coarse_models(self):
-        """Test that MLDA generates the expected exceptions when no coarse_models arg
-        is passed, an empty list is passed or when coarse_models is not a list"""
-        with pytest.raises(TypeError):
-            _, model, _ = mv_simple()
-            with model:
-                MLDA()
-
-        with pytest.raises(ValueError):
-            _, model, _ = mv_simple()
-            with model:
-                MLDA(coarse_models=[])
-
-        with pytest.raises(ValueError):
-            _, model, _ = mv_simple()
-            with model:
-                MLDA(coarse_models=(model, model))
-
-    def test_nonparallelized_chains_are_random(self):
-        """Test that parallel chain are not identical when no parallelisation
-        is applied"""
-        with Model() as coarse_model:
-            Normal("x", 0.3, 1)
-
-        with Model():
-            Normal("x", 0, 1)
-            for stepper in TestMLDA.steppers:
-                step = stepper(coarse_models=[coarse_model])
-                idata = sample(chains=2, cores=1, draws=20, tune=0, step=step, random_seed=1)
-                samples = idata.posterior["x"].values[:, 5]
-                assert len(set(samples)) == 2, f"Non parallelized {stepper} chains are identical."
-
-    def test_parallelized_chains_are_random(self):
-        """Test that parallel chain are
-        not identical when parallelisation
-        is applied"""
-        with Model() as coarse_model:
-            Normal("x", 0.3, 1)
-
-        with Model():
-            Normal("x", 0, 1)
-            for stepper in TestMLDA.steppers:
-                step = stepper(coarse_models=[coarse_model])
-                idata = sample(chains=2, cores=2, draws=20, tune=0, step=step, random_seed=1)
-                samples = idata.posterior["x"].values[:, 5]
-                assert len(set(samples)) == 2, f"Parallelized {stepper} chains are identical."
-
-    def test_acceptance_rate_against_coarseness(self):
-        """Test that the acceptance rate increases
-        when the coarse model is closer to
-        the fine model."""
-        with Model() as coarse_model_0:
-            Normal("x", 5.0, 1.0)
-
-        with Model() as coarse_model_1:
-            Normal("x", 6.0, 2.0)
-
-        with Model() as coarse_model_2:
-            Normal("x", 20.0, 5.0)
-
-        possible_coarse_models = [coarse_model_0, coarse_model_1, coarse_model_2]
-        acc = []
-
-        with Model():
-            Normal("x", 5.0, 1.0)
-            for coarse_model in possible_coarse_models:
-                step = MLDA(coarse_models=[coarse_model], subsampling_rates=3)
-                idata = sample(chains=1, draws=500, tune=100, step=step, random_seed=1)
-                acc.append(idata.sample_stats["accepted"].mean())
-            assert acc[0] > acc[1] > acc[2], (
-                "Acceptance rate is not "
-                "strictly increasing when"
-                "coarse model is closer to "
-                "fine model. Acceptance rates"
-                "were: {}".format(acc)
-            )
-
-    def test_mlda_non_blocked(self):
-        """Test that MLDA correctly creates non-blocked
-        compound steps in level 0 when using a Metropolis
-        base sampler."""
-        _, model = simple_2model_continuous()
-        _, model_coarse = simple_2model_continuous()
-        with model:
-            for stepper in self.steppers:
-                assert isinstance(
-                    stepper(
-                        coarse_models=[model_coarse],
-                        base_sampler="Metropolis",
-                        base_blocked=False,
-                    ).step_method_below,
-                    CompoundStep,
-                )
-
-    def test_mlda_blocked(self):
-        """Test the type of base sampler instantiated
-        when switching base_blocked flag while
-        the base sampler is Metropolis and when
-        the base sampler is DEMetropolisZ."""
-        _, model = simple_2model_continuous()
-        _, model_coarse = simple_2model_continuous()
-        with model:
-            for stepper in self.steppers:
-                assert not isinstance(
-                    stepper(
-                        coarse_models=[model_coarse],
-                        base_sampler="Metropolis",
-                        base_blocked=True,
-                    ).step_method_below,
-                    CompoundStep,
-                )
-                assert isinstance(
-                    stepper(
-                        coarse_models=[model_coarse],
-                        base_sampler="Metropolis",
-                        base_blocked=True,
-                    ).step_method_below,
-                    Metropolis,
-                )
-                assert isinstance(
-                    stepper(coarse_models=[model_coarse]).step_method_below,
-                    DEMetropolisZ,
-                )
-
-    def test_tuning_and_scaling_on(self):
-        """Test that tune and base_scaling change as expected when
-        tuning is on."""
-        np.random.seed(1234)
-        ts = 100
-        _, model = simple_2model_continuous()
-        _, model_coarse = simple_2model_continuous()
-        with model:
-            trace_0 = sample(
-                tune=ts,
-                draws=20,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_sampler="Metropolis",
-                    base_tune_interval=50,
-                    base_scaling=100.0,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=1234,
-                return_inferencedata=False,
-            )
-
-            trace_1 = sample(
-                tune=ts,
-                draws=20,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_tune_target="scaling",
-                    base_tune_interval=50,
-                    base_scaling=100.0,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=1234,
-                return_inferencedata=False,
-            )
-
-            trace_2 = sample(
-                tune=ts,
-                draws=20,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_tune_interval=50,
-                    base_scaling=10,
-                    base_lamb=100.0,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=1234,
-                return_inferencedata=False,
-            )
-
-        assert trace_0.get_sampler_stats("tune", chains=0)[0]
-        assert trace_0.get_sampler_stats("tune", chains=0)[ts - 1]
-        assert not trace_0.get_sampler_stats("tune", chains=0)[ts]
-        assert not trace_0.get_sampler_stats("tune", chains=0)[-1]
-        assert trace_0.get_sampler_stats("base_scaling", chains=0)[0, 0] == 100.0
-        assert trace_0.get_sampler_stats("base_scaling", chains=0)[0, 1] == 100.0
-        assert trace_0.get_sampler_stats("base_scaling", chains=0)[-1, 0] < 100.0
-        assert trace_0.get_sampler_stats("base_scaling", chains=0)[-1, 1] < 100.0
-
-        assert trace_1.get_sampler_stats("tune", chains=0)[0]
-        assert trace_1.get_sampler_stats("tune", chains=0)[ts - 1]
-        assert not trace_1.get_sampler_stats("tune", chains=0)[ts]
-        assert not trace_1.get_sampler_stats("tune", chains=0)[-1]
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[0] == 100.0
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[-1] < 100.0
-
-        assert trace_2.get_sampler_stats("tune", chains=0)[0]
-        assert trace_2.get_sampler_stats("tune", chains=0)[ts - 1]
-        assert not trace_2.get_sampler_stats("tune", chains=0)[ts]
-        assert not trace_2.get_sampler_stats("tune", chains=0)[-1]
-        assert trace_2.get_sampler_stats("base_lambda", chains=0)[0] == 100.0
-        assert trace_2.get_sampler_stats("base_lambda", chains=0)[-1] < 100.0
-
-    def test_tuning_and_scaling_off(self):
-        """Test that tuning is deactivated when sample()'s tune=0 and that
-        MLDA's tune=False is overridden by sample()'s tune."""
-        np.random.seed(12345)
-        _, model = simple_2model_continuous()
-        _, model_coarse = simple_2model_continuous()
-
-        ts_0 = 0
-        with model:
-            trace_0 = sample(
-                tune=ts_0,
-                draws=100,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_sampler="Metropolis",
-                    base_tune_interval=50,
-                    base_scaling=100.0,
-                    tune=False,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=12345,
-                return_inferencedata=False,
-            )
-
-        ts_1 = 100
-        with model:
-            trace_1 = sample(
-                tune=ts_1,
-                draws=20,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_sampler="Metropolis",
-                    base_tune_interval=50,
-                    base_scaling=100.0,
-                    tune=False,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=12345,
-                return_inferencedata=False,
-            )
-
-        assert not trace_0.get_sampler_stats("tune", chains=0)[0]
-        assert not trace_0.get_sampler_stats("tune", chains=0)[-1]
-        assert (
-            trace_0.get_sampler_stats("base_scaling", chains=0)[0, 0]
-            == trace_0.get_sampler_stats("base_scaling", chains=0)[-1, 0]
-            == trace_0.get_sampler_stats("base_scaling", chains=0)[0, 1]
-            == trace_0.get_sampler_stats("base_scaling", chains=0)[-1, 1]
-            == 100.0
-        )
-
-        assert trace_1.get_sampler_stats("tune", chains=0)[0]
-        assert trace_1.get_sampler_stats("tune", chains=0)[ts_1 - 1]
-        assert not trace_1.get_sampler_stats("tune", chains=0)[ts_1]
-        assert not trace_1.get_sampler_stats("tune", chains=0)[-1]
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[0, 0] == 100.0
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[0, 1] == 100.0
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[-1, 0] < 100.0
-        assert trace_1.get_sampler_stats("base_scaling", chains=0)[-1, 1] < 100.0
-
-        ts_2 = 0
-        with model:
-            trace_2 = sample(
-                tune=ts_2,
-                draws=100,
-                step=MLDA(
-                    coarse_models=[model_coarse],
-                    base_tune_interval=50,
-                    base_lamb=100.0,
-                    base_tune_target=None,
-                ),
-                chains=1,
-                discard_tuned_samples=False,
-                random_seed=12345,
-                return_inferencedata=False,
-            )
-
-        assert not trace_2.get_sampler_stats("tune", chains=0)[0]
-        assert not trace_2.get_sampler_stats("tune", chains=0)[-1]
-        assert (
-            trace_2.get_sampler_stats("base_lambda", chains=0)[0]
-            == trace_2.get_sampler_stats("base_lambda", chains=0)[-1]
-            == trace_2.get_sampler_stats("base_lambda", chains=0)[0]
-            == trace_2.get_sampler_stats("base_lambda", chains=0)[-1]
-            == 100.0
-        )
-
-    def test_trace_length(self):
-        """Check if trace length is as expected."""
-        tune = 100
-        draws = 50
-        with Model() as coarse_model:
-            Normal("n", 0, 2.2, size=(3,))
-        with Model():
-            Normal("n", 0, 2, size=(3,))
-            step = MLDA(coarse_models=[coarse_model])
-            idata = sample(tune=tune, draws=draws, step=step, chains=1, discard_tuned_samples=False)
-            assert len(idata.warmup_posterior.draw) == tune
-            assert len(idata.posterior.draw) == draws
-
-    @pytest.mark.parametrize(
-        "variable,has_grad,outcome",
-        [("n", True, 1), ("n", False, 1), ("b", True, 0), ("b", False, 0)],
-    )
-    def test_competence(self, variable, has_grad, outcome):
-        """Test if competence function returns expected
-        results for different models"""
-        with Model() as pmodel:
-            Normal("n", 0, 2, size=(3,))
-            Binomial("b", n=2, p=0.3)
-        assert MLDA.competence(pmodel[variable], has_grad=has_grad) == outcome
-
-    def test_multiple_subsampling_rates(self):
-        """Test that when you give a single integer it is applied to all levels and
-        when you give a list the list is applied correctly."""
-        with Model() as coarse_model_0:
-            Normal("n", 0, 2.2, size=(3,))
-        with Model() as coarse_model_1:
-            Normal("n", 0, 2.1, size=(3,))
-        with Model():
-            Normal("n", 0, 2.0, size=(3,))
-
-            step_1 = MLDA(coarse_models=[coarse_model_0, coarse_model_1], subsampling_rates=3)
-            assert len(step_1.subsampling_rates) == 2
-            assert step_1.subsampling_rates[0] == step_1.subsampling_rates[1] == 3
-
-            step_2 = MLDA(coarse_models=[coarse_model_0, coarse_model_1], subsampling_rates=[3, 4])
-            assert step_2.subsampling_rates[0] == 3
-            assert step_2.subsampling_rates[1] == 4
-
-            with pytest.raises(ValueError):
-                step_3 = MLDA(
-                    coarse_models=[coarse_model_0, coarse_model_1],
-                    subsampling_rates=[3, 4, 10],
-                )
-
-    def test_aem_mu_sigma(self):
-        """Test that AEM estimates mu_B and Sigma_B in
-        the coarse models of a 3-level LR example correctly"""
-        # create data for linear regression
-        if aesara.config.floatX == "float32":
-            p = "float32"
-        else:
-            p = "float64"
-        np.random.seed(123456)
-        size = 200
-        true_intercept = 1
-        true_slope = 2
-        sigma = 1
-        x = np.linspace(0, 1, size, dtype=p)
-        # y = a + b*x
-        true_regression_line = true_intercept + true_slope * x
-        # add noise
-        y = true_regression_line + np.random.normal(0, sigma**2, size)
-        s = np.identity(y.shape[0], dtype=p)
-        np.fill_diagonal(s, sigma**2)
-
-        # forward model Op - here, just the regression equation
-        class ForwardModel(Op):
-            if aesara.config.floatX == "float32":
-                itypes = [at.fvector]
-                otypes = [at.fvector]
-            else:
-                itypes = [at.dvector]
-                otypes = [at.dvector]
-
-            def __init__(self, x, pymc_model):
-                self.x = x
-                self.pymc_model = pymc_model
-
-            def perform(self, node, inputs, outputs):
-                intercept = inputs[0][0]
-                x_coeff = inputs[0][1]
-
-                temp = intercept + x_coeff * x + self.pymc_model.bias.data
-                with self.pymc_model:
-                    set_data({"model_output": temp})
-                outputs[0][0] = np.array(temp)
-
-        # create the coarse models with separate biases
-        mout = []
-        coarse_models = []
-
-        with Model() as coarse_model_0:
-            bias = Data("bias", 3.5 * np.ones(y.shape, dtype=p))
-            mu_B = MutableData("mu_B", -1.3 * np.ones(y.shape, dtype=p))
-            Sigma_B = MutableData("Sigma_B", np.zeros((y.shape[0], y.shape[0]), dtype=p))
-            model_output = MutableData("model_output", np.zeros(y.shape, dtype=p))
-            Sigma_e = Data("Sigma_e", s)
-
-            # Define priors
-            intercept = Normal("Intercept", 0, sigma=20)
-            x_coeff = Normal("x", 0, sigma=20)
-
-            theta = at.as_tensor_variable([intercept, x_coeff])
-
-            mout.append(ForwardModel(x, coarse_model_0))
-
-            # Define likelihood
-            likelihood = MvNormal("y", mu=mout[0](theta) + mu_B, cov=Sigma_e, observed=y)
-
-            coarse_models.append(coarse_model_0)
-
-        with Model() as coarse_model_1:
-            bias = Data("bias", 2.2 * np.ones(y.shape, dtype=p))
-            mu_B = MutableData("mu_B", -2.2 * np.ones(y.shape, dtype=p))
-            Sigma_B = MutableData("Sigma_B", np.zeros((y.shape[0], y.shape[0]), dtype=p))
-            model_output = MutableData("model_output", np.zeros(y.shape, dtype=p))
-            Sigma_e = Data("Sigma_e", s)
-
-            # Define priors
-            intercept = Normal("Intercept", 0, sigma=20)
-            x_coeff = Normal("x", 0, sigma=20)
-
-            theta = at.as_tensor_variable([intercept, x_coeff])
-
-            mout.append(ForwardModel(x, coarse_model_1))
-
-            # Define likelihood
-            likelihood = MvNormal("y", mu=mout[1](theta) + mu_B, cov=Sigma_e, observed=y)
-
-            coarse_models.append(coarse_model_1)
-
-        # fine model and inference
-        with Model() as model:
-            bias = Data("bias", np.zeros(y.shape, dtype=p))
-            model_output = MutableData("model_output", np.zeros(y.shape, dtype=p))
-            Sigma_e = Data("Sigma_e", s)
-
-            # Define priors
-            intercept = Normal("Intercept", 0, sigma=20)
-            x_coeff = Normal("x", 0, sigma=20)
-
-            theta = at.as_tensor_variable([intercept, x_coeff])
-
-            mout.append(ForwardModel(x, model))
-
-            # Define likelihood
-            likelihood = MvNormal("y", mu=mout[-1](theta), cov=Sigma_e, observed=y)
-
-            step_mlda = MLDA(coarse_models=coarse_models, adaptive_error_model=True)
-
-            trace_mlda = sample(
-                draws=100,
-                step=step_mlda,
-                chains=1,
-                tune=200,
-                discard_tuned_samples=True,
-                random_seed=84759238,
-            )
-
-            m0 = step_mlda.step_method_below.model_below.mu_B.get_value()
-            s0 = step_mlda.step_method_below.model_below.Sigma_B.get_value()
-            m1 = step_mlda.model_below.mu_B.get_value()
-            s1 = step_mlda.model_below.Sigma_B.get_value()
-
-            assert np.allclose(m0, -3.5)
-            assert np.allclose(m1, -2.2)
-            assert np.allclose(s0, 0, atol=1e-3)
-            assert np.allclose(s1, 0, atol=1e-3)
-
-    def test_variance_reduction(self):
-        """
-        Test if the right stats are outputed when variance reduction is used in MLDA,
-        if the output estimates are close (VR estimate vs. standard estimate from
-        the first chain) and if the variance of VR is lower. Uses a linear regression
-        model with multiple levels where approximate levels have fewer data.
-
-        """
-        # arithmetic precision
-        if aesara.config.floatX == "float32":
-            p = "float32"
-        else:
-            p = "float64"
-
-        # set up the model and data
-        seed = 12345
-        np.random.seed(seed)
-        size = 100
-        true_intercept = 1
-        true_slope = 2
-        sigma = 0.1
-        x = np.linspace(0, 1, size, dtype=p)
-        # y = a + b*x
-        true_regression_line = true_intercept + true_slope * x
-        # add noise
-        y = true_regression_line + np.random.normal(0, sigma**2, size)
-        s = sigma
-
-        x_coarse_0 = x[::3]
-        y_coarse_0 = y[::3]
-        x_coarse_1 = x[::2]
-        y_coarse_1 = y[::2]
-
-        # MCMC parameters
-        ndraws = 200
-        ntune = 100
-        nsub = 3
-        nchains = 1
-
-        # define likelihoods with different Q
-        class Likelihood(Op):
-            if aesara.config.floatX == "float32":
-                itypes = [at.fvector]
-                otypes = [at.fscalar]
-            else:
-                itypes = [at.dvector]
-                otypes = [at.dscalar]
-
-            def __init__(self, x, y, pymc_model):
-                self.x = x
-                self.y = y
-                self.pymc_model = pymc_model
-
-            def perform(self, node, inputs, outputs):
-                intercept = inputs[0][0]
-                x_coeff = inputs[0][1]
-
-                temp = np.array(intercept + x_coeff * self.x, dtype=p)
-                with self.pymc_model:
-                    set_data({"Q": np.array(x_coeff, dtype=p)})
-                outputs[0][0] = np.array(
-                    -(0.5 / s**2) * np.sum((temp - self.y) ** 2, dtype=p), dtype=p
-                )
-
-        # run four MLDA steppers for all combinations of
-        # base_sampler and forward model
-        for stepper in ["Metropolis", "DEMetropolisZ"]:
-            mout = []
-            coarse_models = []
-
-            with Model() as coarse_model_0:
-                if aesara.config.floatX == "float32":
-                    Q = MutableData("Q", np.float32(0.0))
-                else:
-                    Q = MutableData("Q", np.float64(0.0))
-
-                # Define priors
-                intercept = Normal("Intercept", true_intercept, sigma=1)
-                x_coeff = Normal("x", true_slope, sigma=1)
-
-                theta = at.as_tensor_variable([intercept, x_coeff])
-
-                mout.append(Likelihood(x_coarse_0, y_coarse_0, coarse_model_0))
-                Potential("likelihood", mout[0](theta))
-
-                coarse_models.append(coarse_model_0)
-
-            with Model() as coarse_model_1:
-                if aesara.config.floatX == "float32":
-                    Q = MutableData("Q", np.float32(0.0))
-                else:
-                    Q = MutableData("Q", np.float64(0.0))
-
-                # Define priors
-                intercept = Normal("Intercept", true_intercept, sigma=1)
-                x_coeff = Normal("x", true_slope, sigma=1)
-
-                theta = at.as_tensor_variable([intercept, x_coeff])
-
-                mout.append(Likelihood(x_coarse_1, y_coarse_1, coarse_model_1))
-                Potential("likelihood", mout[1](theta))
-
-                coarse_models.append(coarse_model_1)
-
-            with Model() as model:
-                if aesara.config.floatX == "float32":
-                    Q = MutableData("Q", np.float32(0.0))
-                else:
-                    Q = MutableData("Q", np.float64(0.0))
-
-                # Define priors
-                intercept = Normal("Intercept", true_intercept, sigma=1)
-                x_coeff = Normal("x", true_slope, sigma=1)
-
-                theta = at.as_tensor_variable([intercept, x_coeff])
-
-                mout.append(Likelihood(x, y, model))
-                Potential("likelihood", mout[-1](theta))
-
-                step = MLDA(
-                    coarse_models=coarse_models,
-                    base_sampler=stepper,
-                    subsampling_rates=nsub,
-                    variance_reduction=True,
-                    store_Q_fine=True,
-                )
-
-                trace = sample(
-                    draws=ndraws,
-                    step=step,
-                    chains=nchains,
-                    tune=ntune,
-                    cores=1,
-                    discard_tuned_samples=True,
-                    random_seed=seed,
-                    return_inferencedata=False,
-                )
-
-                # get fine level stats (standard method)
-                Q_2 = trace.get_sampler_stats("Q_2").reshape((nchains, ndraws))
-                Q_mean_standard = Q_2.mean(axis=1).mean()
-                Q_se_standard = np.sqrt(Q_2.var() / az.ess(np.array(Q_2, np.float64)))
-
-                # get VR stats
-                Q_mean_vr, Q_se_vr = extract_Q_estimate(trace, 3)
-
-                # check that returned values are floats and finite.
-                assert isinstance(Q_mean_standard, np.floating)
-                assert np.isfinite(Q_mean_standard)
-                assert isinstance(Q_mean_vr, np.floating)
-                assert np.isfinite(Q_mean_vr)
-                assert isinstance(Q_se_standard, np.floating)
-                assert np.isfinite(Q_se_standard)
-                assert isinstance(Q_se_vr, np.floating)
-                assert np.isfinite(Q_se_vr)
-
-                # check consistency of QoI across levels.
-                Q_1_0 = np.concatenate(trace.get_sampler_stats("Q_1_0")).reshape(
-                    (nchains, ndraws * nsub)
-                )
-                Q_2_1 = np.concatenate(trace.get_sampler_stats("Q_2_1")).reshape((nchains, ndraws))
-                # This used to be a scrict zero equality!
-                assert np.isclose(Q_1_0.mean(axis=1), 0.0, atol=1e-4)
-                assert np.isclose(Q_2_1.mean(axis=1), 0.0, atol=1e-4)
-
-
 class TestRVsAssignmentSteps:
     """
     Test that step methods convert input RVs to respective value vars
@@ -1378,7 +661,6 @@ class TestRVsAssignmentSteps:
             (Slice, {}),
             (DEMetropolis, {}),
             (DEMetropolisZ, {}),
-            # (MLDA, {}),  # TODO
         ],
     )
     def test_continuous_steps(self, step, step_kwargs):
diff --git a/pymc/tests/test_types.py b/pymc/tests/test_types.py
index b08f0d26b..f5e884abf 100644
--- a/pymc/tests/test_types.py
+++ b/pymc/tests/test_types.py
@@ -20,7 +20,7 @@ import numpy as np
 from pymc.distributions import Normal
 from pymc.model import Model
 from pymc.sampling import sample
-from pymc.step_methods import MLDA, NUTS, HamiltonianMC, Metropolis, Slice
+from pymc.step_methods import NUTS, HamiltonianMC, Metropolis, Slice
 
 
 class TestType:
@@ -59,39 +59,3 @@ class TestType:
         for sampler in self.samplers:
             with model:
                 sample(draws=10, tune=10, chains=1, step=sampler())
-
-    @aesara.config.change_flags({"floatX": "float64", "warn_float64": "ignore"})
-    def test_float64_MLDA(self):
-        data = np.random.randn(5)
-
-        with Model() as coarse_model:
-            x = Normal("x", initval=np.array(1.0, dtype="float64"))
-            obs = Normal("obs", mu=x, sigma=1.0, observed=data + 0.5)
-
-        with Model() as model:
-            x = Normal("x", initval=np.array(1.0, dtype="float64"))
-            obs = Normal("obs", mu=x, sigma=1.0, observed=data)
-
-        assert x.dtype == "float64"
-        assert obs.dtype == "float64"
-
-        with model:
-            sample(draws=10, tune=10, chains=1, step=MLDA(coarse_models=[coarse_model]))
-
-    @aesara.config.change_flags({"floatX": "float32", "warn_float64": "warn"})
-    def test_float32_MLDA(self):
-        data = np.random.randn(5).astype("float32")
-
-        with Model() as coarse_model:
-            x = Normal("x", initval=np.array(1.0, dtype="float32"))
-            obs = Normal("obs", mu=x, sigma=1.0, observed=data + 0.5)
-
-        with Model() as model:
-            x = Normal("x", initval=np.array(1.0, dtype="float32"))
-            obs = Normal("obs", mu=x, sigma=1.0, observed=data)
-
-        assert x.dtype == "float32"
-        assert obs.dtype == "float32"
-
-        with model:
-            sample(draws=10, tune=10, chains=1, step=MLDA(coarse_models=[coarse_model]))
diff --git a/pymc/tests/test_util.py b/pymc/tests/test_util.py
index 345469a6f..1a0cfe1d1 100644
--- a/pymc/tests/test_util.py
+++ b/pymc/tests/test_util.py
@@ -14,13 +14,20 @@
 
 import numpy as np
 import pytest
+import xarray
 
 from cachetools import cached
 
 import pymc as pm
 
 from pymc.distributions.transforms import RVTransform
-from pymc.util import UNSET, hash_key, hashable, locally_cachedmethod
+from pymc.util import (
+    UNSET,
+    dataset_to_point_list,
+    hash_key,
+    hashable,
+    locally_cachedmethod,
+)
 
 
 class TestTransformName:
@@ -142,3 +149,18 @@ def test_unset_repr(capsys):
     help(fn)
     captured = capsys.readouterr()
     assert "a=UNSET" in captured.out
+
+
+def test_dataset_to_point_list():
+    ds = xarray.Dataset()
+    ds["A"] = xarray.DataArray([[1, 2, 3]] * 2, dims=("chain", "draw"))
+    pl = dataset_to_point_list(ds)
+    assert isinstance(pl, list)
+    assert len(pl) == 6
+    assert isinstance(pl[0], dict)
+    assert isinstance(pl[0]["A"], np.ndarray)
+
+    # Check that non-str keys are caught
+    ds[3] = xarray.DataArray([1, 2, 3])
+    with pytest.raises(ValueError, match="must be str"):
+        dataset_to_point_list(ds)