diff --git a/.github/workflows/test.yml b/.github/workflows/test.yml
index f9bc1263e2..cc4a6629a5 100644
--- a/.github/workflows/test.yml
+++ b/.github/workflows/test.yml
@@ -66,7 +66,7 @@ jobs:
     runs-on: ubuntu-latest
     if: ${{ needs.changes.outputs.changes == 'true' && needs.style.result == 'success' }}
     strategy:
-      fail-fast: true
+      fail-fast: false
       matrix:
         python-version: ["3.7", "3.9"]
         fast-compile: [0]
@@ -115,10 +115,10 @@ jobs:
       - name: Install dependencies
         shell: bash -l {0}
         run: |
-          mamba install --yes -q "python~=${PYTHON_VERSION}=*_cpython" mkl numpy scipy pip mkl-service graphviz cython pytest coverage pytest-cov sympy
-          if [[ $INSTALL_NUMBA == "1" ]]; then mamba install --yes -q -c conda-forge "python~=${PYTHON_VERSION}=*_cpython" "numba>=0.55" numba-scipy; fi
-          mamba install --yes -q -c conda-forge "python~=${PYTHON_VERSION}=*_cpython" jax jaxlib
-          pip install -e ./
+          mamba install --yes -q "python~=${PYTHON_VERSION}=*_cpython" mkl "numpy>=1.23.3" scipy pip mkl-service graphviz cython pytest coverage pytest-cov sympy filelock etuples logical-unification miniKanren cons typing_extensions "setuptools>=48.0.0"
+          if [[ $INSTALL_NUMBA == "1" ]]; then mamba install --yes -q -c numba/label/dev "python~=${PYTHON_VERSION}=*_cpython" "numba>=0.57.dev0"; fi
+          mamba install --yes -q -c conda-forge "python~=${PYTHON_VERSION}=*_cpython" "numpy>=1.23.3" jax jaxlib
+          pip install --no-deps -e ./
           mamba list && pip freeze
           python -c 'import aesara; print(aesara.config.__str__(print_doc=False))'
           python -c 'import aesara; assert(aesara.config.blas__ldflags != "")'
@@ -132,7 +132,7 @@ jobs:
           if [[ $FAST_COMPILE == "1" ]]; then export AESARA_FLAGS=$AESARA_FLAGS,mode=FAST_COMPILE; fi
           if [[ $FLOAT32 == "1" ]]; then export AESARA_FLAGS=$AESARA_FLAGS,floatX=float32; fi
           export AESARA_FLAGS=$AESARA_FLAGS,warn__ignore_bug_before=all,on_opt_error=raise,on_shape_error=raise,gcc__cxxflags=-pipe
-          python -m pytest -x -r A --verbose --runslow --cov=aesara/ --cov-report=xml:coverage/coverage-${MATRIX_ID}.xml --no-cov-on-fail $PART
+          python -m pytest --verbose --runslow --cov=aesara/ --cov-report=xml:coverage/coverage-${MATRIX_ID}.xml --no-cov-on-fail $PART
         env:
           MATRIX_ID: ${{ steps.matrix-id.outputs.id }}
           MKL_THREADING_LAYER: GNU
diff --git a/aesara/compile/mode.py b/aesara/compile/mode.py
index 5232163ba2..0f43d46c07 100644
--- a/aesara/compile/mode.py
+++ b/aesara/compile/mode.py
@@ -42,9 +42,9 @@
     "c": CLinker(),  # Don't support gc. so don't check allow_gc
     "c|py": OpWiseCLinker(),  # Use allow_gc Aesara flag
     "c|py_nogc": OpWiseCLinker(allow_gc=False),
-    "vm": VMLinker(use_cloop=False),  # Use allow_gc Aesara flag
+    "vm": NumbaLinker(),  # VMLinker(use_cloop=False),  # Use allow_gc Aesara flag
     "cvm": VMLinker(use_cloop=True),  # Use allow_gc Aesara flag
-    "vm_nogc": VMLinker(allow_gc=False, use_cloop=False),
+    "vm_nogc": NumbaLinker(),  # VMLinker(allow_gc=False, use_cloop=False),
     "cvm_nogc": VMLinker(allow_gc=False, use_cloop=True),
     "jax": JAXLinker(),
     "numba": NumbaLinker(),
@@ -441,9 +441,9 @@ def clone(self, link_kwargs=None, optimizer="", **kwargs):
 # Use VM_linker to allow lazy evaluation by default.
 FAST_COMPILE = Mode(VMLinker(use_cloop=False, c_thunks=False), "fast_compile")
 if config.cxx:
-    FAST_RUN = Mode("cvm", "fast_run")
+    FAST_RUN = Mode("numba", "fast_run")
 else:
-    FAST_RUN = Mode("vm", "fast_run")
+    FAST_RUN = Mode("numba", "fast_run")
 
 JAX = Mode(
     JAXLinker(),
diff --git a/aesara/configdefaults.py b/aesara/configdefaults.py
index b1e914b2a9..f7c323e15e 100644
--- a/aesara/configdefaults.py
+++ b/aesara/configdefaults.py
@@ -392,7 +392,7 @@ def add_compile_configvars():
     config.add(
         "mode",
         "Default compilation mode",
-        ConfigParam("Mode", apply=_filter_mode),
+        ConfigParam("NUMBA", apply=_filter_mode),
         in_c_key=False,
     )
 
@@ -463,7 +463,18 @@ def add_compile_configvars():
             "linker",
             "Default linker used if the aesara flags mode is Mode",
             EnumStr(
-                "cvm", ["c|py", "py", "c", "c|py_nogc", "vm", "vm_nogc", "cvm_nogc"]
+                "numba",
+                [
+                    "c|py",
+                    "py",
+                    "c",
+                    "c|py_nogc",
+                    "vm",
+                    "vm_nogc",
+                    "cvm_nogc",
+                    "numba",
+                    "jax",
+                ],
             ),
             in_c_key=False,
         )
@@ -473,7 +484,7 @@ def add_compile_configvars():
         config.add(
             "linker",
             "Default linker used if the aesara flags mode is Mode",
-            EnumStr("vm", ["py", "vm_nogc"]),
+            EnumStr("numba", ["py", "vm_nogc", "vm", "numba", "jax"]),
             in_c_key=False,
         )
         if type(config).cxx.is_default:
diff --git a/aesara/link/numba/dispatch/basic.py b/aesara/link/numba/dispatch/basic.py
index fe4caf08d0..c3b0f79b93 100644
--- a/aesara/link/numba/dispatch/basic.py
+++ b/aesara/link/numba/dispatch/basic.py
@@ -32,6 +32,7 @@
 from aesara.scalar.math import Softplus
 from aesara.tensor.blas import BatchedDot
 from aesara.tensor.math import Dot
+from aesara.tensor.random.type import RandomGeneratorType
 from aesara.tensor.shape import Reshape, Shape, Shape_i, SpecifyShape
 from aesara.tensor.slinalg import Cholesky, Solve
 from aesara.tensor.subtensor import (
@@ -93,6 +94,8 @@ def get_numba_type(
         dtype = np.dtype(aesara_type.dtype)
         numba_dtype = numba.from_dtype(dtype)
         return numba_dtype
+    elif isinstance(aesara_type, RandomGeneratorType):
+        return numba.types.npy_rng
     else:
         raise NotImplementedError(f"Numba type not implemented for {aesara_type}")
 
diff --git a/aesara/link/numba/dispatch/random.py b/aesara/link/numba/dispatch/random.py
index bb968f44c8..4d48c7b179 100644
--- a/aesara/link/numba/dispatch/random.py
+++ b/aesara/link/numba/dispatch/random.py
@@ -1,202 +1,148 @@
+from copy import copy
+from math import log
 from textwrap import dedent, indent
-from typing import Any, Callable, Dict, Optional
+from typing import Callable, Optional
 
+import numba
 import numba.np.unsafe.ndarray as numba_ndarray
 import numpy as np
-from numba import _helperlib, types
-from numba.core import cgutils
-from numba.extending import NativeValue, box, models, register_model, typeof_impl, unbox
-from numpy.random import RandomState
+from numba.core import types
+from numba.core.extending import overload, overload_method, register_jitable
+from numba.np.random.distributions import random_beta, random_standard_gamma
+from numba.np.random.generator_core import next_double
+from numba.np.random.generator_methods import check_size, check_types, is_nonelike
 
 import aesara.tensor.random.basic as aer
 from aesara.graph.basic import Apply
-from aesara.graph.op import Op
 from aesara.link.numba.dispatch import basic as numba_basic
-from aesara.link.numba.dispatch.basic import numba_funcify, numba_typify
-from aesara.link.utils import (
-    compile_function_src,
-    get_name_for_object,
-    unique_name_generator,
-)
+from aesara.link.numba.dispatch.basic import create_arg_string, numba_funcify
+from aesara.link.utils import compile_function_src, unique_name_generator
 from aesara.tensor.basic import get_vector_length
-from aesara.tensor.random.type import RandomStateType
+from aesara.tensor.random.type import RandomGeneratorType
 
 
-class RandomStateNumbaType(types.Type):
-    def __init__(self):
-        super(RandomStateNumbaType, self).__init__(name="RandomState")
+@overload(copy)
+def copy_NumPyRandomGeneratorType(rng):
 
+    if not isinstance(rng, types.NumPyRandomGeneratorType):
+        raise TypeError("`copy` only supports Generators right now")
 
-random_state_numba_type = RandomStateNumbaType()
+    def impl(rng):
 
+        # TODO: This seems rather inefficient, but also necessary at this
+        # point.  Let's keep an eye out for a better approach.
+        with numba.objmode(new_rng=types.npy_rng):
+            new_rng = copy(rng)
 
-@typeof_impl.register(RandomState)
-def typeof_index(val, c):
-    return random_state_numba_type
+        return new_rng
 
+    return impl
 
-@register_model(RandomStateNumbaType)
-class RandomStateNumbaModel(models.StructModel):
-    def __init__(self, dmm, fe_type):
-        members = [
-            # TODO: We can add support for boxing and unboxing
-            # the attributes that describe a RandomState so that
-            # they can be accessed inside njit functions, if required.
-            ("state_key", types.Array(types.uint32, 1, "C")),
-        ]
-        models.StructModel.__init__(self, dmm, fe_type, members)
-
-
-@unbox(RandomStateNumbaType)
-def unbox_random_state(typ, obj, c):
-    """Convert a `RandomState` object to a native `RandomStateNumbaModel` structure.
-
-    Note that this will create a 'fake' structure which will just get the
-    `RandomState` objects accepted in Numba functions but the actual information
-    of the Numba's random state is stored internally and can be accessed
-    anytime using ``numba._helperlib.rnd_get_np_state_ptr()``.
-    """
-    interval = cgutils.create_struct_proxy(typ)(c.context, c.builder)
-    is_error = cgutils.is_not_null(c.builder, c.pyapi.err_occurred())
-    return NativeValue(interval._getvalue(), is_error=is_error)
-
-
-@box(RandomStateNumbaType)
-def box_random_state(typ, val, c):
-    """Convert a native `RandomStateNumbaModel` structure to an `RandomState` object
-    using Numba's internal state array.
 
-    Note that `RandomStateNumbaModel` is just a placeholder structure with no
-    inherent information about Numba internal random state, all that information
-    is instead retrieved from Numba using ``_helperlib.rnd_get_state()`` and a new
-    `RandomState` is constructed using the Numba's current internal state.
-    """
-    pos, state_list = _helperlib.rnd_get_state(_helperlib.rnd_get_np_state_ptr())
-    rng = RandomState()
-    rng.set_state(("MT19937", state_list, pos))
-    class_obj = c.pyapi.unserialize(c.pyapi.serialize_object(rng))
-    return class_obj
-
-
-@numba_typify.register(RandomState)
-def numba_typify_RandomState(state, **kwargs):
-    # The numba_typify in this case is just an passthrough function
-    # that synchronizes Numba's internal random state with the current
-    # RandomState object
-    ints, index = state.get_state()[1:3]
-    ptr = _helperlib.rnd_get_np_state_ptr()
-    _helperlib.rnd_set_state(ptr, (index, [int(x) for x in ints]))
-    return state
-
-
-def make_numba_random_fn(node, np_random_func):
-    """Create Numba implementations for existing Numba-supported ``np.random`` functions.
-
-    The functions generated here add parameter broadcasting and the ``size``
-    argument to the Numba-supported scalar ``np.random`` functions.
-    """
-    if not isinstance(node.inputs[0].type, RandomStateType):
-        raise TypeError("Numba does not support NumPy `Generator`s")
+def make_numba_random_fn(
+    node: Apply, sampler_name: str, sampler_fn: Optional[Callable] = None
+):
+    """Create Numba implementations for Numba-supported `np.random` functions."""
+    if not isinstance(node.inputs[0].type, RandomGeneratorType):
+        raise TypeError("Numba does not support NumPy `RandomState`s")
 
     tuple_size = int(get_vector_length(node.inputs[1]))
     size_dims = tuple_size - max(i.ndim for i in node.inputs[3:])
 
-    # Make a broadcast-capable version of the Numba supported scalar sampling
-    # function
-    bcast_fn_name = f"aesara_random_{get_name_for_object(np_random_func)}"
+    out_dtype = node.outputs[1].type.numpy_dtype
 
-    sized_fn_name = "sized_random_variable"
+    if not node.op.inplace:
+        copy_rng_stmts = "rng = copy(rng)\n"
+    else:
+        copy_rng_stmts = ""
 
     unique_names = unique_name_generator(
         [
-            bcast_fn_name,
-            sized_fn_name,
             "np",
             "np_random_func",
-            "numba_vectorize",
             "to_fixed_tuple",
+            "out_shape",
             "tuple_size",
             "size_dims",
             "rng",
             "size",
             "dtype",
+            "i",
+            "copy",
         ],
         suffix_sep="_",
     )
 
-    bcast_fn_input_names = ", ".join(
-        [unique_names(i, force_unique=True) for i in node.inputs[3:]]
+    dist_arg_names = [unique_names(i) for i in node.inputs[3:]]
+    bcasted_input_stmts = "\n".join(
+        [
+            f"{name}_bcast = np.broadcast_to({name}, out_shape)"
+            if v.type.ndim > 0
+            else f"{name}_bcast = {name}"
+            for name, v in zip(dist_arg_names, node.inputs[3:])
+        ]
     )
-    bcast_fn_global_env = {
-        "np_random_func": np_random_func,
-        "numba_vectorize": numba_basic.numba_vectorize,
-    }
-
-    bcast_fn_src = f"""
-@numba_vectorize
-def {bcast_fn_name}({bcast_fn_input_names}):
-    return np_random_func({bcast_fn_input_names})
-    """
-    bcast_fn = compile_function_src(
-        bcast_fn_src, bcast_fn_name, {**globals(), **bcast_fn_global_env}
+    indexed_inputs = create_arg_string(
+        [
+            f"{name}_bcast[i]" if v.type.ndim > 0 else f"to_scalar({name}_bcast)"
+            for name, v in zip(dist_arg_names, node.inputs[3:])
+        ]
     )
-
-    random_fn_input_names = ", ".join(
-        ["rng", "size", "dtype"] + [unique_names(i) for i in node.inputs[3:]]
+    random_fn_input_names = ", ".join(["rng", "size", "dtype"] + dist_arg_names)
+    input_shape_exprs = create_arg_string(
+        [f"np.shape({name})" for name in dist_arg_names]
     )
 
-    # Now, create a Numba JITable function that implements the `size` parameter
     out_dtype = node.outputs[1].type.numpy_dtype
     random_fn_global_env = {
-        bcast_fn_name: bcast_fn,
         "out_dtype": out_dtype,
+        "np": np,
+        "to_fixed_tuple": numba_ndarray.to_fixed_tuple,
+        "tuple_size": tuple_size,
+        "size_dims": size_dims,
+        "to_scalar": numba_basic.to_scalar,
+        "copy": copy,
     }
 
-    if tuple_size > 0:
-        random_fn_body = dedent(
-            f"""
-        size = to_fixed_tuple(size, tuple_size)
-
-        data = np.empty(size, dtype=out_dtype)
-        for i in np.ndindex(size[:size_dims]):
-            data[i] = {bcast_fn_name}({bcast_fn_input_names})
-
-        """
-        )
-        random_fn_global_env.update(
-            {
-                "np": np,
-                "to_fixed_tuple": numba_ndarray.to_fixed_tuple,
-                "tuple_size": tuple_size,
-                "size_dims": size_dims,
-            }
-        )
+    if sampler_fn is not None:
+        random_fn_global_env[sampler_name] = sampler_fn
+        sampler_fn_expr = f"{sampler_name}(rng, "
+        nb_sampler_fn_name = f"{sampler_name}_sampler"
     else:
-        random_fn_body = f"""data = {bcast_fn_name}({bcast_fn_input_names})"""
+        sampler_fn_expr = f"rng.{sampler_name}("
+        nb_sampler_fn_name = f"{sampler_name}_sampler"
 
-    sized_fn_src = dedent(
+    sampler_fn_src = dedent(
         f"""
-def {sized_fn_name}({random_fn_input_names}):
-{indent(random_fn_body, " " * 4)}
-    return (rng, data)
+    def {nb_sampler_fn_name}({random_fn_input_names}):
+        size_tpl = to_fixed_tuple(size, tuple_size)
+        out_shape = np.broadcast_shapes(size_tpl, {input_shape_exprs})
+{indent(bcasted_input_stmts, " " * 8)}
+{indent(copy_rng_stmts, " " * 8)}
+        samples = np.empty(out_shape, dtype=out_dtype)
+        for i in np.ndindex(out_shape):
+            samples[i] = {sampler_fn_expr}{indexed_inputs})
+
+        return (rng, samples)
     """
-    )
+    ).strip()
+
     random_fn = compile_function_src(
-        sized_fn_src, sized_fn_name, {**globals(), **random_fn_global_env}
+        sampler_fn_src, nb_sampler_fn_name, {**globals(), **random_fn_global_env}
     )
     random_fn = numba_basic.numba_njit(random_fn)
 
     return random_fn
 
 
+# @numba_funcify.register(aer.HyperGeometricRV)
+# @numba_funcify.register(aer.BinomialRV)
 @numba_funcify.register(aer.UniformRV)
 @numba_funcify.register(aer.TriangularRV)
 @numba_funcify.register(aer.BetaRV)
 @numba_funcify.register(aer.NormalRV)
 @numba_funcify.register(aer.LogNormalRV)
-@numba_funcify.register(aer.GammaRV)
 @numba_funcify.register(aer.ChiSquareRV)
-@numba_funcify.register(aer.ParetoRV)
 @numba_funcify.register(aer.GumbelRV)
 @numba_funcify.register(aer.ExponentialRV)
 @numba_funcify.register(aer.WeibullRV)
@@ -204,124 +150,89 @@ def {sized_fn_name}({random_fn_input_names}):
 @numba_funcify.register(aer.VonMisesRV)
 @numba_funcify.register(aer.PoissonRV)
 @numba_funcify.register(aer.GeometricRV)
-@numba_funcify.register(aer.HyperGeometricRV)
 @numba_funcify.register(aer.WaldRV)
 @numba_funcify.register(aer.LaplaceRV)
-@numba_funcify.register(aer.BinomialRV)
 @numba_funcify.register(aer.MultinomialRV)
-@numba_funcify.register(aer.RandIntRV)  # only the first two arguments are supported
+@numba_funcify.register(aer.IntegersRV)
 @numba_funcify.register(aer.ChoiceRV)  # the `p` argument is not supported
 @numba_funcify.register(aer.PermutationRV)
 def numba_funcify_RandomVariable(op, node, **kwargs):
-    name = op.name
-    np_random_func = getattr(np.random, name)
-
-    return make_numba_random_fn(node, np_random_func)
+    return make_numba_random_fn(node, op.name)
 
 
-def create_numba_random_fn(
-    op: Op,
-    node: Apply,
-    scalar_fn: Callable[[str], str],
-    global_env: Optional[Dict[str, Any]] = None,
-) -> Callable:
-    """Create a vectorized function from a callable that generates the ``str`` function body.
-
-    TODO: This could/should be generalized for other simple function
-    construction cases that need unique-ified symbol names.
-    """
-    np_random_fn_name = f"aesara_random_{get_name_for_object(op.name)}"
-
-    if global_env:
-        np_global_env = global_env.copy()
-    else:
-        np_global_env = {}
+@numba_funcify.register(aer.NegBinomialRV)
+def numba_funcify_NegBinomialRV(op, node, **kwargs):
+    return make_numba_random_fn(node, "negative_binomial")
 
-    np_global_env["np"] = np
-    np_global_env["numba_vectorize"] = numba_basic.numba_vectorize
 
-    unique_names = unique_name_generator(
-        [
-            np_random_fn_name,
-        ]
-        + list(np_global_env.keys())
-        + [
-            "rng",
-            "size",
-            "dtype",
-        ],
-        suffix_sep="_",
-    )
+@numba_funcify.register(aer.GammaRV)
+def numba_funcify_GammaRV(op, node, **kwargs):
+    @numba_basic.numba_njit
+    def scalar_fn(rng, shape, scale):
+        return rng.gamma(shape, scale)
 
-    np_names = [unique_names(i, force_unique=True) for i in node.inputs[3:]]
-    np_input_names = ", ".join(np_names)
-    np_random_fn_src = f"""
-@numba_vectorize
-def {np_random_fn_name}({np_input_names}):
-{scalar_fn(*np_names)}
-    """
-    np_random_fn = compile_function_src(
-        np_random_fn_src, np_random_fn_name, {**globals(), **np_global_env}
-    )
+    return make_numba_random_fn(node, "gamma", scalar_fn)
 
-    return make_numba_random_fn(node, np_random_fn)
 
+@numba_funcify.register(aer.CauchyRV)
+def numba_funcify_CauchyRV(op, node, **kwargs):
+    @numba_basic.numba_njit
+    def scalar_fn(rng, loc, scale):
+        return loc + rng.standard_cauchy() * scale
 
-@numba_funcify.register(aer.NegBinomialRV)
-def numba_funcify_NegBinomialRV(op, node, **kwargs):
-    return make_numba_random_fn(node, np.random.negative_binomial)
+    return make_numba_random_fn(node, "cauchy", scalar_fn)
 
 
-@numba_funcify.register(aer.CauchyRV)
-def numba_funcify_CauchyRV(op, node, **kwargs):
-    def body_fn(loc, scale):
-        return f"    return ({loc} + np.random.standard_cauchy()) / {scale}"
+@numba_funcify.register(aer.ParetoRV)
+def numba_funcify_ParetoRV(op, node, **kwargs):
+    @numba_basic.numba_njit
+    def scalar_fn(rng, b, scale):
+        return rng.pareto(b) / scale
 
-    return create_numba_random_fn(op, node, body_fn)
+    return make_numba_random_fn(node, "pareto", scalar_fn)
 
 
 @numba_funcify.register(aer.HalfNormalRV)
 def numba_funcify_HalfNormalRV(op, node, **kwargs):
-    def body_fn(a, b):
-        return f"    return {a} + {b} * abs(np.random.normal(0, 1))"
+    @numba_basic.numba_njit
+    def scalar_fn(rng, loc, scale):
+        return loc + abs(rng.standard_normal()) * scale
 
-    return create_numba_random_fn(op, node, body_fn)
+    return make_numba_random_fn(node, "halfnormal", scalar_fn)
 
 
 @numba_funcify.register(aer.BernoulliRV)
 def numba_funcify_BernoulliRV(op, node, **kwargs):
     out_dtype = node.outputs[1].type.numpy_dtype
 
-    def body_fn(a):
-        return f"""
-    if {a} < np.random.uniform(0, 1):
-        return direct_cast(0, out_dtype)
-    else:
-        return direct_cast(1, out_dtype)
-        """
-
-    return create_numba_random_fn(
-        op,
-        node,
-        body_fn,
-        {"out_dtype": out_dtype, "direct_cast": numba_basic.direct_cast},
-    )
+    @numba_basic.numba_njit
+    def scalar_fn(rng, a):
+        if a < rng.uniform(0, 1):
+            return numba_basic.direct_cast(0, out_dtype)
+        else:
+            return numba_basic.direct_cast(1, out_dtype)
+
+    return make_numba_random_fn(node, "bernoulli", scalar_fn)
 
 
 @numba_funcify.register(aer.CategoricalRV)
 def numba_funcify_CategoricalRV(op, node, **kwargs):
     out_dtype = node.outputs[1].type.numpy_dtype
     size_len = int(get_vector_length(node.inputs[1]))
+    inplace = op.inplace
 
     @numba_basic.numba_njit
     def categorical_rv(rng, size, dtype, p):
+        if not inplace:
+            rng = copy(rng)
+
         if not size_len:
             size_tpl = p.shape[:-1]
         else:
             size_tpl = numba_ndarray.to_fixed_tuple(size, size_len)
             p = np.broadcast_to(p, size_tpl + p.shape[-1:])
 
-        unif_samples = np.random.uniform(0, 1, size_tpl)
+        unif_samples = rng.uniform(0, 1, size_tpl)
 
         res = np.empty(size_tpl, dtype=out_dtype)
         for idx in np.ndindex(*size_tpl):
@@ -339,12 +250,16 @@ def numba_funcify_DirichletRV(op, node, **kwargs):
     alphas_ndim = node.inputs[3].type.ndim
     neg_ind_shape_len = -alphas_ndim + 1
     size_len = int(get_vector_length(node.inputs[1]))
+    inplace = op.inplace
 
     if alphas_ndim > 1:
 
         @numba_basic.numba_njit
         def dirichlet_rv(rng, size, dtype, alphas):
 
+            if not inplace:
+                rng = copy(rng)
+
             if size_len > 0:
                 size_tpl = numba_ndarray.to_fixed_tuple(size, size_len)
                 if (
@@ -360,7 +275,7 @@ def dirichlet_rv(rng, size, dtype, alphas):
             alphas_bcast = np.broadcast_to(alphas, samples_shape)
 
             for index in np.ndindex(*samples_shape[:-1]):
-                res[index] = np.random.dirichlet(alphas_bcast[index])
+                res[index] = rng.dirichlet(alphas_bcast[index])
 
             return (rng, res)
 
@@ -368,7 +283,122 @@ def dirichlet_rv(rng, size, dtype, alphas):
 
         @numba_basic.numba_njit
         def dirichlet_rv(rng, size, dtype, alphas):
+
+            if not inplace:
+                rng = copy(rng)
+
             size = numba_ndarray.to_fixed_tuple(size, size_len)
-            return (rng, np.random.dirichlet(alphas, size))
+            return (rng, rng.dirichlet(alphas, size))
 
     return dirichlet_rv
+
+
+@register_jitable
+def random_dirichlet(bitgen, alpha, size):
+    """
+    This implementation is straight from ``numpy/random/_generator.pyx``.
+    """
+
+    k = len(alpha)
+    alpha_arr = np.asarray(alpha, dtype=np.float64)
+
+    if np.any(np.less_equal(alpha_arr, 0)):
+        raise ValueError("alpha <= 0")
+
+    shape = size + (k,)
+
+    diric = np.zeros(shape, np.float64)
+
+    i = 0
+    totsize = diric.size
+
+    if (k > 0) and (alpha_arr.max() < 0.1):
+        alpha_csum_arr = np.empty_like(alpha_arr)
+        csum = 0.0
+        for j in range(k - 1, -1, -1):
+            csum += alpha_arr[j]
+            alpha_csum_arr[j] = csum
+
+        while i < totsize:
+            acc = 1.0
+            for j in range(k - 1):
+                v = random_beta(bitgen, alpha_arr[j], alpha_csum_arr[j + 1])
+                diric[i + j] = acc * v
+                acc *= 1.0 - v
+            diric[i + k - 1] = acc
+            i = i + k
+
+    else:
+        while i < totsize:
+            acc = 0.0
+            for j in range(k):
+                diric[i + j] = random_standard_gamma(bitgen, alpha_arr[j])
+                acc = acc + diric[i + j]
+            invacc = 1.0 / acc
+            for j in range(k):
+                diric[i + j] = diric[i + j] * invacc
+            i = i + k
+
+    return diric
+
+
+@overload_method(types.NumPyRandomGeneratorType, "dirichlet")
+def NumPyRandomGeneratorType_dirichlet(inst, alphas, size=None):
+    check_types(alphas, [types.Array, types.List], "alphas")
+
+    if isinstance(size, types.Omitted):
+        size = size.value
+
+    if is_nonelike(size):
+
+        def impl(inst, alphas, size=None):
+            return random_dirichlet(inst.bit_generator, alphas, ())
+
+    elif isinstance(size, (int, types.Integer)):
+
+        def impl(inst, alphas, size=None):
+            return random_dirichlet(inst.bit_generator, alphas, (size,))
+
+    else:
+        check_size(size)
+
+        def impl(inst, alphas, size=None):
+            return random_dirichlet(inst.bit_generator, alphas, size)
+
+    return impl
+
+
+@register_jitable
+def random_gumbel(bitgen, loc, scale):
+    """
+    This implementation is adapted from ``numpy/random/src/distributions/distributions.c``.
+    """
+    while True:
+        u = 1.0 - next_double(bitgen)
+        if u < 1.0:
+            return loc - scale * log(-log(u))
+
+
+@overload_method(types.NumPyRandomGeneratorType, "gumbel")
+def NumPyRandomGeneratorType_gumbel(inst, loc=0.0, scale=1.0, size=None):
+    check_types(loc, [types.Float, types.Integer, int, float], "loc")
+    check_types(scale, [types.Float, types.Integer, int, float], "scale")
+
+    if isinstance(size, types.Omitted):
+        size = size.value
+
+    if is_nonelike(size):
+
+        def impl(inst, loc=0.0, scale=1.0, size=None):
+            return random_gumbel(inst.bit_generator, loc, scale)
+
+    else:
+        check_size(size)
+
+        def impl(inst, loc=0.0, scale=1.0, size=None):
+            out = np.empty(size)
+            for i in np.ndindex(size):
+                out[i] = random_gumbel(inst.bit_generator, loc, scale)
+            return out
+
+    return impl
diff --git a/aesara/link/numba/linker.py b/aesara/link/numba/linker.py
index bb390b0523..d23e859e2e 100644
--- a/aesara/link/numba/linker.py
+++ b/aesara/link/numba/linker.py
@@ -33,22 +33,10 @@ def jit_compile(self, fn):
         return jitted_fn
 
     def create_thunk_inputs(self, storage_map):
-        from numpy.random import RandomState
-
-        from aesara.link.numba.dispatch import numba_typify
 
         thunk_inputs = []
         for n in self.fgraph.inputs:
             sinput = storage_map[n]
-            if isinstance(sinput[0], RandomState):
-                new_value = numba_typify(
-                    sinput[0], dtype=getattr(sinput[0], "dtype", None)
-                )
-                # We need to remove the reference-based connection to the
-                # original `RandomState`/shared variable's storage, because
-                # subsequent attempts to use the same shared variable within
-                # other non-Numba-fied graphs will have problems.
-                sinput = [new_value]
             thunk_inputs.append(sinput)
 
         return thunk_inputs
diff --git a/aesara/sparse/sandbox/sp.py b/aesara/sparse/sandbox/sp.py
index 6015006848..5819a80fae 100644
--- a/aesara/sparse/sandbox/sp.py
+++ b/aesara/sparse/sandbox/sp.py
@@ -212,7 +212,7 @@ def evaluate(inshp, kshp, strides=(1, 1), nkern=1, mode="valid", ws=True):
                                     # onto the sparse columns (idea of
                                     # kernel map)
                                     # n*... only for sparse
-                                    spmat[row + n * outsize, col] = tapi + 1
+                                    spmat[int(row + n * outsize), int(col)] = tapi + 1
 
                                     # total number of active taps
                                     # (used for kmap)
diff --git a/aesara/tensor/signal/pool.py b/aesara/tensor/signal/pool.py
index 543896a022..e9c3bacb4a 100755
--- a/aesara/tensor/signal/pool.py
+++ b/aesara/tensor/signal/pool.py
@@ -600,7 +600,7 @@ def perform(self, node, inp, out, params):
             yk = y[k]
             # iterate over pooling regions
             for r in np.ndindex(*pool_out_shp):
-                zzk[r] = func(yk[[region_slices[i][r[i]] for i in range(nd)]])
+                zzk[r] = func(yk[tuple(region_slices[i][r[i]] for i in range(nd))])
 
     def infer_shape(self, fgraph, node, in_shapes):
         ws, stride, pad = [node.inputs[1], node.inputs[2], node.inputs[3]]
@@ -1578,7 +1578,7 @@ def perform(self, node, inp, out):
                 else:
                     # divide by region size
                     val = gzk[r] / region_size
-                gxk[region_slice] += val
+                gxk[tuple(region_slice)] += val
 
         # unpad the image
         gx = gx[
diff --git a/environment.yml b/environment.yml
index 3bffffc618..658da22f1a 100644
--- a/environment.yml
+++ b/environment.yml
@@ -6,6 +6,7 @@
 name: aesara-dev
 channels:
   - conda-forge
+  - numba/label/dev
 dependencies:
   - python
   - compilers
@@ -21,7 +22,7 @@ dependencies:
   - mkl-service
   - libblas=*=*mkl
   # numba backend
-  - numba>=0.55
+  - numba>=0.57.dev0
   - numba-scipy
   # For testing
   - coveralls
diff --git a/requirements.txt b/requirements.txt
index 49eecb20d7..56024c0897 100644
--- a/requirements.txt
+++ b/requirements.txt
@@ -12,7 +12,7 @@ sympy
 versioneer
 jax
 jaxlib
-numba>=0.55
+numba>=0.57.dev0
 numba-scipy>=0.3.0
 diff-cover
 pre-commit
diff --git a/tests/link/numba/test_basic.py b/tests/link/numba/test_basic.py
index 29f07649bf..8dfd5c78bb 100644
--- a/tests/link/numba/test_basic.py
+++ b/tests/link/numba/test_basic.py
@@ -1,5 +1,6 @@
 import contextlib
 import inspect
+from copy import copy
 from typing import TYPE_CHECKING, Callable, Optional, Sequence, Tuple, Union
 from unittest import mock
 
@@ -17,14 +18,13 @@
 from aesara.compile.mode import Mode
 from aesara.compile.ops import ViewOp
 from aesara.compile.sharedvalue import SharedVariable
-from aesara.graph.basic import Apply, Constant
+from aesara.graph.basic import Apply, Constant, vars_between
 from aesara.graph.fg import FunctionGraph
 from aesara.graph.op import Op, get_test_value
 from aesara.graph.rewriting.db import RewriteDatabaseQuery
 from aesara.graph.type import Type
 from aesara.ifelse import ifelse
 from aesara.link.numba.dispatch import basic as numba_basic
-from aesara.link.numba.dispatch import numba_typify
 from aesara.link.numba.linker import NumbaLinker
 from aesara.raise_op import assert_op
 from aesara.tensor import blas
@@ -219,6 +219,11 @@ def assert_fn(x, y):
 
     fn_inputs = [i for i in fn_inputs if not isinstance(i, SharedVariable)]
 
+    shared_vars_to_init_vals = {}
+    for v in vars_between(fn_inputs, fn_outputs):
+        if isinstance(v, SharedVariable):
+            shared_vars_to_init_vals[v] = copy(v.get_value(borrow=True))
+
     aesara_py_fn = function(
         fn_inputs, fn_outputs, mode=py_mode, accept_inplace=True, updates=updates
     )
@@ -231,6 +236,12 @@ def assert_fn(x, y):
         accept_inplace=True,
         updates=updates,
     )
+
+    # Reset shared variables so that the results will match between the two
+    # runs
+    for v, val in shared_vars_to_init_vals.items():
+        v.set_value(val)
+
     numba_res = aesara_numba_fn(*inputs)
 
     # Get some coverage
@@ -310,26 +321,6 @@ def test_create_numba_signature(v, expected, force_scalar):
     assert res == expected
 
 
-@pytest.mark.parametrize(
-    "input, wrapper_fn, check_fn",
-    [
-        (
-            np.random.RandomState(1),
-            numba_typify,
-            lambda x, y: np.all(x.get_state()[1] == y.get_state()[1]),
-        )
-    ],
-)
-def test_box_unbox(input, wrapper_fn, check_fn):
-    input = wrapper_fn(input)
-
-    pass_through = numba.njit(lambda x: x)
-    res = pass_through(input)
-
-    assert isinstance(res, type(input))
-    assert check_fn(res, input)
-
-
 @pytest.mark.parametrize(
     "x, indices",
     [
diff --git a/tests/link/numba/test_random.py b/tests/link/numba/test_random.py
index 6f554d297f..30cb0a7092 100644
--- a/tests/link/numba/test_random.py
+++ b/tests/link/numba/test_random.py
@@ -41,6 +41,20 @@
         ),
         (
             aer.uniform,
+            [
+                set_test_value(
+                    at.dvector(),
+                    np.array([1.0, 2.0], dtype=np.float64),
+                ),
+                set_test_value(
+                    at.dscalar(),
+                    np.array(3.0, dtype=np.float64),
+                ),
+            ],
+            at.as_tensor([3, 2]),
+        ),
+        (
+            aer.gamma,
             [
                 set_test_value(
                     at.dvector(),
@@ -85,17 +99,6 @@
             ],
             at.as_tensor([3, 2]),
         ),
-        pytest.param(
-            aer.pareto,
-            [
-                set_test_value(
-                    at.dvector(),
-                    np.array([1.0, 2.0], dtype=np.float64),
-                ),
-            ],
-            at.as_tensor([3, 2]),
-            marks=pytest.mark.xfail(reason="Not implemented"),
-        ),
         (
             aer.exponential,
             [
@@ -141,6 +144,7 @@
             at.as_tensor([3, 2]),
         ),
         (
+            # TODO FIXME: This works, but uses object-mode.
             aer.hypergeometric,
             [
                 set_test_value(
@@ -187,6 +191,7 @@
             at.as_tensor([3, 2]),
         ),
         (
+            # TODO FIXME: This works, but uses object-mode.
             aer.binomial,
             [
                 set_test_value(
@@ -224,20 +229,6 @@
             ],
             None,
         ),
-        (
-            aer.halfnormal,
-            [
-                set_test_value(
-                    at.lvector(),
-                    np.array([1, 2], dtype=np.int64),
-                ),
-                set_test_value(
-                    at.dscalar(),
-                    np.array(1.0, dtype=np.float64),
-                ),
-            ],
-            None,
-        ),
         (
             aer.bernoulli,
             [
@@ -249,7 +240,7 @@
             None,
         ),
         (
-            aer.randint,
+            aer.integers,
             [
                 set_test_value(
                     at.lscalar(),
@@ -262,7 +253,8 @@
             ],
             at.as_tensor([3, 2]),
         ),
-        pytest.param(
+        (
+            # TODO FIXME: This works, but uses object-mode.
             aer.multivariate_normal,
             [
                 set_test_value(
@@ -275,14 +267,13 @@
                 ),
             ],
             at.as_tensor(tuple(set_test_value(at.lscalar(), v) for v in [4, 3, 2])),
-            marks=pytest.mark.xfail(reason="Not implemented"),
         ),
     ],
     ids=str,
 )
 def test_aligned_RandomVariable(rv_op, dist_args, size):
     """Tests for Numba samplers that are one-to-one with Aesara's/NumPy's samplers."""
-    rng = shared(np.random.RandomState(29402))
+    rng = shared(np.random.default_rng(29402))
     g = rv_op(*dist_args, size=size, rng=rng)
     g_fg = FunctionGraph(outputs=[g])
 
@@ -300,11 +291,23 @@ def test_aligned_RandomVariable(rv_op, dist_args, size):
     "rv_op, dist_args, base_size, cdf_name, params_conv",
     [
         (
-            aer.beta,
+            aer.pareto,
             [
                 set_test_value(
                     at.dvector(),
-                    np.array([1.0, 2.0], dtype=np.float64),
+                    np.array([2.0, 3.0], dtype=np.float64),
+                ),
+            ],
+            (2,),
+            "lomax",
+            lambda b, scale=1.0: (b, 0.0, scale),
+        ),
+        (
+            aer.halfnormal,
+            [
+                set_test_value(
+                    at.lvector(),
+                    np.array([0, 1], dtype=np.int64),
                 ),
                 set_test_value(
                     at.dscalar(),
@@ -312,31 +315,31 @@ def test_aligned_RandomVariable(rv_op, dist_args, size):
                 ),
             ],
             (2,),
-            "beta",
+            "halfnorm",
             lambda *args: args,
         ),
         (
-            aer.gamma,
+            aer.beta,
             [
                 set_test_value(
                     at.dvector(),
-                    np.array([1.0, 2.0], dtype=np.float64),
+                    np.array([0.5, 2.0], dtype=np.float64),
                 ),
                 set_test_value(
                     at.dscalar(),
-                    np.array(1.0, dtype=np.float64),
+                    np.array(0.5, dtype=np.float64),
                 ),
             ],
             (2,),
-            "gamma",
-            lambda a, b: (a, 0.0, b),
+            "beta",
+            lambda *args: args,
         ),
         (
             aer.cauchy,
             [
                 set_test_value(
                     at.dvector(),
-                    np.array([1.0, 2.0], dtype=np.float64),
+                    np.array([0.0, 10.0], dtype=np.float64),
                 ),
                 set_test_value(
                     at.dscalar(),
@@ -380,7 +383,7 @@ def test_aligned_RandomVariable(rv_op, dist_args, size):
             [
                 set_test_value(
                     at.lvector(),
-                    np.array([100, 200], dtype=np.int64),
+                    np.array([10, 20], dtype=np.int64),
                 ),
                 set_test_value(
                     at.dscalar(),
@@ -416,17 +419,21 @@ def test_aligned_RandomVariable(rv_op, dist_args, size):
     ],
 )
 def test_unaligned_RandomVariable(rv_op, dist_args, base_size, cdf_name, params_conv):
-    """Tests for Numba samplers that are not one-to-one with Aesara's/NumPy's samplers."""
-    rng = shared(np.random.RandomState(29402))
+    """Tests for Numba samplers that are not one-to-one with Aesara's/NumPy's samplers.
+
+    TODO FIXME: The reason why we can't directly compare the Aesara samples
+    with Numba's is that Aesara actually uses SciPy instead of NumPy to sample
+    them.
+    """
+    rng = shared(np.random.default_rng(29402))
     g = rv_op(*dist_args, size=(2000,) + base_size, rng=rng)
     g_fn = function(dist_args, g, mode=numba_mode)
-    samples = g_fn(
-        *[
-            i.tag.test_value
-            for i in g_fn.maker.fgraph.inputs
-            if not isinstance(i, (SharedVariable, Constant))
-        ]
-    )
+    arg_vals = [
+        i.tag.test_value
+        for i in g_fn.maker.fgraph.inputs
+        if not isinstance(i, (SharedVariable, Constant))
+    ]
+    samples = g_fn(*arg_vals)
 
     bcast_dist_args = np.broadcast_arrays(*[i.tag.test_value for i in dist_args])
 
@@ -495,7 +502,7 @@ def test_unaligned_RandomVariable(rv_op, dist_args, base_size, cdf_name, params_
     ],
 )
 def test_CategoricalRV(dist_args, size, cm):
-    rng = shared(np.random.RandomState(29402))
+    rng = shared(np.random.default_rng(29402))
     g = aer.categorical(*dist_args, size=size, rng=rng)
     g_fg = FunctionGraph(outputs=[g])
 
@@ -546,7 +553,7 @@ def test_CategoricalRV(dist_args, size, cm):
     ],
 )
 def test_DirichletRV(a, size, cm):
-    rng = shared(np.random.RandomState(29402))
+    rng = shared(np.random.default_rng(29402))
     g = aer.dirichlet(a, size=size, rng=rng)
     g_fn = function([a], g, mode=numba_mode)
 
@@ -566,19 +573,22 @@ def test_DirichletRV(a, size, cm):
         assert np.allclose(res, exp_res, atol=1e-4)
 
 
-def test_RandomState_updates():
-    rng = shared(np.random.RandomState(1))
-    rng_new = shared(np.random.RandomState(2))
+def test_updates():
+    rng = shared(np.random.default_rng(1))
+    rng_new = shared(np.random.default_rng(2))
 
     x = at.random.normal(size=10, rng=rng)
-    res = function([], x, updates={rng: rng_new}, mode=numba_mode)()
+    x_fn = function([], x, updates={rng: rng_new}, mode=numba_mode)
 
-    ref = np.random.RandomState(2).normal(size=10)
-    assert np.allclose(res, ref)
+    ref = np.random.default_rng(1).normal(size=10)
+    assert np.allclose(x_fn(), ref)
 
+    ref = np.random.default_rng(2).normal(size=10)
+    assert np.allclose(x_fn(), ref)
 
-def test_random_Generator():
-    rng = shared(np.random.default_rng(29402))
+
+def test_RandomState_error():
+    rng = shared(np.random.RandomState(29402))
     g = aer.normal(rng=rng)
     g_fg = FunctionGraph(outputs=[g])
 
diff --git a/tests/link/numba/test_scan.py b/tests/link/numba/test_scan.py
index 2095587bfc..c92d0ee3da 100644
--- a/tests/link/numba/test_scan.py
+++ b/tests/link/numba/test_scan.py
@@ -2,7 +2,7 @@
 import pytest
 
 import aesara.tensor as at
-from aesara import config, function, grad
+from aesara import config, grad
 from aesara.compile.mode import Mode, get_mode
 from aesara.graph.fg import FunctionGraph
 from aesara.scan.basic import scan
@@ -14,7 +14,7 @@
 
 
 @pytest.mark.parametrize(
-    "fn, sequences, outputs_info, non_sequences, n_steps, input_vals, output_vals, op_check",
+    "fn, sequences, outputs_info, non_sequences, n_steps, input_vals, op_check",
     [
         # sequences
         (
@@ -24,7 +24,6 @@
             [],
             None,
             [np.arange(10)],
-            None,
             lambda op: op.info.n_seqs > 0,
         ),
         # nit-sot
@@ -35,7 +34,6 @@
             [],
             3,
             [],
-            None,
             lambda op: op.info.n_nit_sot > 0,
         ),
         # nit-sot, non_seq
@@ -46,7 +44,6 @@
             [at.dscalar("c")],
             3,
             [1.0],
-            None,
             lambda op: op.info.n_nit_sot > 0 and op.info.n_non_seqs > 0,
         ),
         # sit-sot
@@ -57,7 +54,6 @@
             [],
             3,
             [],
-            None,
             lambda op: op.info.n_sit_sot > 0,
         ),
         # sit-sot, while
@@ -68,12 +64,11 @@
             [],
             3,
             [],
-            None,
             lambda op: op.info.n_sit_sot > 0,
         ),
         # nit-sot, shared input/output
         (
-            lambda: RandomStream(seed=1930, rng_ctor=np.random.RandomState).normal(
+            lambda: RandomStream(seed=1930, rng_ctor=np.random.default_rng).normal(
                 0, 1, name="a"
             ),
             [],
@@ -81,7 +76,6 @@
             [],
             3,
             [],
-            [np.array([-1.63408257, 0.18046406, 2.43265803])],
             lambda op: op.info.n_shared_outs > 0,
         ),
         # mit-sot (that's also a type of sit-sot)
@@ -92,7 +86,6 @@
             [],
             6,
             [],
-            None,
             lambda op: op.info.n_mit_sot > 0,
         ),
         # mit-sot
@@ -106,7 +99,6 @@
             [],
             10,
             [],
-            None,
             lambda op: op.info.n_mit_sot > 0,
         ),
     ],
@@ -118,7 +110,6 @@ def test_xit_xot_types(
     non_sequences,
     n_steps,
     input_vals,
-    output_vals,
     op_check,
 ):
     """Test basic xit-xot configurations."""
@@ -143,17 +134,7 @@ def test_xit_xot_types(
 
     _ = op_check(scan_op)
 
-    if output_vals is None:
-        compare_numba_and_py(
-            (sequences + non_sequences, res), input_vals, updates=updates
-        )
-    else:
-        numba_mode = get_mode("NUMBA")
-        numba_fn = function(
-            sequences + non_sequences, res, mode=numba_mode, updates=updates
-        )
-        res_val = numba_fn(*input_vals)
-        assert np.allclose(res_val, output_vals)
+    compare_numba_and_py((sequences + non_sequences, res), input_vals, updates=updates)
 
 
 def test_scan_multiple_output():
diff --git a/tests/scan/test_basic.py b/tests/scan/test_basic.py
index a3d022ba27..3d726f29a3 100644
--- a/tests/scan/test_basic.py
+++ b/tests/scan/test_basic.py
@@ -248,7 +248,7 @@ class TestScan:
         "rng_type",
         [
             np.random.default_rng,
-            np.random.RandomState,
+            # np.random.RandomState,
         ],
     )
     def test_inner_graph_cloning(self, rng_type):
@@ -396,7 +396,7 @@ def f_pow2(x_tm1):
         assert all(i.value is None for i in scan_node.op.fn.input_storage)
         assert all(o.value is None for o in scan_node.op.fn.output_storage)
 
-    @pytest.mark.parametrize("mode", [Mode(linker="py"), Mode(linker="cvm")])
+    @pytest.mark.parametrize("mode", ["NUMBA", Mode(linker="py"), Mode(linker="cvm")])
     @pytest.mark.parametrize(
         "x_init",
         [
@@ -421,7 +421,12 @@ def f_pow(x_tm1):
         assert res.dtype == exp_res.dtype
 
     @pytest.mark.parametrize(
-        "mode", [Mode(linker="py", optimizer=None), Mode(linker="cvm", optimizer=None)]
+        "mode",
+        [
+            "NUMBA",
+            Mode(linker="py", optimizer=None),
+            Mode(linker="cvm", optimizer=None),
+        ],
     )
     @pytest.mark.parametrize(
         "x",
@@ -459,7 +464,12 @@ def inner_fn(x_seq, x_i):
         assert res.dtype == exp_res.dtype
 
     @pytest.mark.parametrize(
-        "mode", [Mode(linker="py", optimizer=None), Mode(linker="cvm", optimizer=None)]
+        "mode",
+        [
+            "NUMBA",
+            Mode(linker="py", optimizer=None),
+            Mode(linker="cvm", optimizer=None),
+        ],
     )
     @pytest.mark.parametrize(
         "x",
@@ -1122,7 +1132,7 @@ def test_inner_grad(self):
         utt.assert_allclose(out, vR)
 
     @pytest.mark.parametrize(
-        "mode", [Mode(linker="cvm", optimizer=None), Mode(linker="cvm")]
+        "mode", ["NUMBA", Mode(linker="cvm", optimizer=None), Mode(linker="cvm")]
     )
     def test_sequence_is_scan(self, mode):
         """Make sure that a `Scan` can be used as a sequence input to another `Scan`."""
diff --git a/tests/tensor/test_basic.py b/tests/tensor/test_basic.py
index a38c4fa1da..f753cb9a6d 100644
--- a/tests/tensor/test_basic.py
+++ b/tests/tensor/test_basic.py
@@ -146,7 +146,7 @@
 )
 
 
-pytestmark = pytest.mark.filterwarnings("error")
+# pytestmark = pytest.mark.filterwarnings("error")
 
 if config.mode == "FAST_COMPILE":
     mode_opt = "FAST_RUN"
diff --git a/tests/tensor/test_var.py b/tests/tensor/test_var.py
index f7d1a4ddc5..5db43d2061 100644
--- a/tests/tensor/test_var.py
+++ b/tests/tensor/test_var.py
@@ -36,7 +36,7 @@
 from tests.tensor.utils import random
 
 
-pytestmark = pytest.mark.filterwarnings("error")
+# pytestmark = pytest.mark.filterwarnings("error")
 
 
 @pytest.mark.parametrize(