diff --git a/docs/manuals/measurement/readout_mitigation.rst b/docs/manuals/measurement/readout_mitigation.rst
index c0672a6118..80e74f65fa 100644
--- a/docs/manuals/measurement/readout_mitigation.rst
+++ b/docs/manuals/measurement/readout_mitigation.rst
@@ -47,25 +47,16 @@ experiments to generate the corresponding mitigators.
     import numpy as np
     import matplotlib.pyplot as plt
     from qiskit import QuantumCircuit
+    from qiskit.quantum_info import Operator
     from qiskit.visualization import plot_distribution
+    from qiskit_experiments.data_processing import LocalReadoutMitigator
     from qiskit_experiments.library import LocalReadoutError, CorrelatedReadoutError
 
     from qiskit_aer import AerSimulator
     from qiskit_ibm_runtime.fake_provider import FakePerth
 
-    from qiskit.result.mitigation.utils import (
-        expval_with_stddev,
-        str2diag,
-        counts_probability_vector
-    )
-
     backend = AerSimulator.from_backend(FakePerth())
 
-.. jupyter-execute::
-
-    shots = 1024
-    qubits = [0,1,2,3]
-    num_qubits = len(qubits)
 
 Standard mitigation experiment
 ------------------------------
@@ -76,13 +67,14 @@ circuits, one for all “0” and one for all “1” results.
 
 .. jupyter-execute::
 
+    shots = 1024
+    qubits = [0,1,2,3]
+    num_qubits = len(qubits)
+
     exp = LocalReadoutError(qubits)
     for c in exp.circuits():
         print(c)
 
-
-.. jupyter-execute::
-
     exp.analysis.set_options(plot=True)
     result = exp.run(backend)
     mitigator = result.analysis_results("Local Readout Mitigator").value
@@ -102,9 +94,9 @@ The individual mitigation matrices can be read off the mitigator.
 
 .. jupyter-execute::
 
-    for m in mitigator._mitigation_mats:
-        print(m)
-        print()
+    for qubit in mitigator.qubits:
+        print(f"Qubit: {qubit}")
+        print(mitigator.mitigation_matrix(qubits=qubit))
 
 
 Mitigation example
@@ -118,13 +110,9 @@ Mitigation example
         qc.cx(i - 1, i)
     qc.measure_all()
 
-.. jupyter-execute::
-
     counts = backend.run(qc, shots=shots, seed_simulator=42, method="density_matrix").result().get_counts()
     unmitigated_probs = {label: count / shots for label, count in counts.items()}
 
-.. jupyter-execute::
-
     mitigated_quasi_probs = mitigator.quasi_probabilities(counts)
     mitigated_stddev = mitigated_quasi_probs._stddev_upper_bound
     mitigated_probs = (mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities())
@@ -144,15 +132,20 @@ Expectation value
 .. jupyter-execute::
 
     diagonal_labels = ["ZZZZ", "ZIZI", "IZII", "1ZZ0"]
-    ideal_expectation = []
-    diagonals = [str2diag(d) for d in diagonal_labels]
+    diagonals = [
+        np.diag(np.real(Operator.from_label(d).to_matrix()))
+        for d in diagonal_labels
+    ]
+
+    # Create a mitigator with no mitigation so that we can use its
+    # expectation_values method to generate an unmitigated expectation value to
+    # compare to the mitigated one.
+    identity_mitigator = LocalReadoutMitigator([np.eye(2) for _ in range(4)])
+
     qubit_index = {i: i for i in range(num_qubits)}
-    unmitigated_probs_vector, _ = counts_probability_vector(unmitigated_probs, qubit_index=qubit_index)
-    unmitigated_expectation = [expval_with_stddev(d, unmitigated_probs_vector, shots) for d in diagonals]
+    unmitigated_expectation = [identity_mitigator.expectation_value(counts, d) for d in diagonals]
     mitigated_expectation = [mitigator.expectation_value(counts, d) for d in diagonals]
 
-.. jupyter-execute::
-
     mitigated_expectation_values, mitigated_stddev = zip(*mitigated_expectation)
     unmitigated_expectation_values, unmitigated_stddev = zip(*unmitigated_expectation)
     legend = ['Mitigated Expectation', 'Unmitigated Expectation']
diff --git a/qiskit_experiments/data_processing/__init__.py b/qiskit_experiments/data_processing/__init__.py
index 56a05ddfa3..b075c1effe 100644
--- a/qiskit_experiments/data_processing/__init__.py
+++ b/qiskit_experiments/data_processing/__init__.py
@@ -82,6 +82,15 @@
     BaseDiscriminator
     SkLDA
     SkQDA
+
+Mitigators
+==========
+.. autosummary::
+    :toctree: ../stubs/
+
+    BaseReadoutMitigator
+    LocalReadoutMitigator
+    CorrelatedReadoutMitigator
 """
 
 from .data_action import DataAction, TrainableDataAction
@@ -104,4 +113,7 @@
 
 from .data_processor import DataProcessor
 from .discriminator import BaseDiscriminator
+from .mitigation.base_readout_mitigator import BaseReadoutMitigator
+from .mitigation.correlated_readout_mitigator import CorrelatedReadoutMitigator
+from .mitigation.local_readout_mitigator import LocalReadoutMitigator
 from .sklearn_discriminators import SkLDA, SkQDA
diff --git a/qiskit_experiments/data_processing/mitigation/__init__.py b/qiskit_experiments/data_processing/mitigation/__init__.py
new file mode 100644
index 0000000000..5cdc4ebf1d
--- /dev/null
+++ b/qiskit_experiments/data_processing/mitigation/__init__.py
@@ -0,0 +1,22 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2021.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+
+"""Readout error mitigation."""
+from .base_readout_mitigator import BaseReadoutMitigator
+from .correlated_readout_mitigator import CorrelatedReadoutMitigator
+from .local_readout_mitigator import LocalReadoutMitigator
+from .utils import (
+    counts_probability_vector,
+    expval_with_stddev,
+    stddev,
+    str2diag,
+)
diff --git a/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py
new file mode 100644
index 0000000000..c37fa333ea
--- /dev/null
+++ b/qiskit_experiments/data_processing/mitigation/base_readout_mitigator.py
@@ -0,0 +1,81 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Base class for readout error mitigation.
+"""
+
+from abc import ABC, abstractmethod
+from typing import Optional, List, Iterable, Tuple, Union, Callable
+
+import numpy as np
+
+from qiskit.result.counts import Counts
+from qiskit.result.distributions.quasi import QuasiDistribution
+
+
+class BaseReadoutMitigator(ABC):
+    """Base readout error mitigator class."""
+
+    @abstractmethod
+    def quasi_probabilities(
+        self,
+        data: Counts,
+        qubits: Iterable[int] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> QuasiDistribution:
+        """Convert counts to a dictionary of quasi-probabilities
+
+        Args:
+            data: Counts to be mitigated.
+            qubits: the physical qubits measured to obtain the counts clbits.
+                If None these are assumed to be qubits [0, ..., N-1]
+                for N-bit counts.
+            clbits: Optional, marginalize counts to just these bits.
+            shots: Optional, the total number of shots, if None shots will
+                be calculated as the sum of all counts.
+
+        Returns:
+            QuasiDistribution: A dictionary containing pairs of [output, mean] where "output"
+                is the key in the dictionaries,
+                which is the length-N bitstring of a measured standard basis state,
+                and "mean" is the mean of non-zero quasi-probability estimates.
+        """
+
+    @abstractmethod
+    def expectation_value(
+        self,
+        data: Counts,
+        diagonal: Union[Callable, dict, str, np.ndarray],
+        qubits: Iterable[int] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> Tuple[float, float]:
+        """Calculate the expectation value of a diagonal Hermitian operator.
+
+        Args:
+            data: Counts object to be mitigated.
+            diagonal: the diagonal operator. This may either be specified
+                      as a string containing I,Z,0,1 characters, or as a
+                      real valued 1D array_like object supplying the full diagonal,
+                      or as a dictionary, or as Callable.
+            qubits: the physical qubits measured to obtain the counts clbits.
+                    If None these are assumed to be qubits [0, ..., N-1]
+                    for N-bit counts.
+            clbits: Optional, marginalize counts to just these bits.
+            shots: Optional, the total number of shots, if None shots will
+                be calculated as the sum of all counts.
+
+        Returns:
+            The mean and an upper bound of the standard deviation of operator
+            expectation value calculated from the current counts.
+        """
diff --git a/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py
new file mode 100644
index 0000000000..8099347285
--- /dev/null
+++ b/qiskit_experiments/data_processing/mitigation/correlated_readout_mitigator.py
@@ -0,0 +1,270 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Readout mitigator class based on the A-matrix inversion method
+"""
+
+import math
+from typing import Optional, List, Tuple, Iterable, Callable, Union, Dict
+import numpy as np
+
+from qiskit.exceptions import QiskitError
+
+from qiskit.result.distributions.quasi import QuasiDistribution
+from qiskit.result.counts import Counts
+from .base_readout_mitigator import BaseReadoutMitigator
+from .utils import counts_probability_vector, z_diagonal, str2diag
+
+
+class CorrelatedReadoutMitigator(BaseReadoutMitigator):
+    """N-qubit readout error mitigator.
+
+    Mitigates :meth:`expectation_value` and :meth:`quasi_probabilities`.
+    The mitigation_matrix should be calibrated using qiskit experiments.
+    This mitigation method should be used in case the readout errors of the qubits
+    are assumed to be correlated. The mitigation_matrix of *N* qubits is of size
+    :math:`2^N x 2^N` so the mitigation complexity is :math:`O(4^N)`.
+    """
+
+    def __init__(self, assignment_matrix: np.ndarray, qubits: Optional[Iterable[int]] = None):
+        """Initialize a CorrelatedReadoutMitigator
+
+        Args:
+            assignment_matrix: readout error assignment matrix.
+            qubits: Optional, the measured physical qubits for mitigation.
+
+        Raises:
+            QiskitError: matrix size does not agree with number of qubits
+        """
+        if np.any(assignment_matrix < 0) or not np.allclose(np.sum(assignment_matrix, axis=0), 1):
+            raise QiskitError("Assignment matrix columns must be valid probability distributions")
+        assignment_matrix = np.asarray(assignment_matrix, dtype=float)
+        matrix_qubits_num = int(math.log2(assignment_matrix.shape[0]))
+        if qubits is None:
+            self._num_qubits = matrix_qubits_num
+            self._qubits = range(self._num_qubits)
+        else:
+            if len(qubits) != matrix_qubits_num:
+                raise QiskitError(
+                    f"The number of given qubits ({len(qubits)}) is different than the number of "
+                    f"qubits inferred from the matrices ({matrix_qubits_num})"
+                )
+            self._qubits = qubits
+            self._num_qubits = len(self._qubits)
+        self._qubit_index = dict(zip(self._qubits, range(self._num_qubits)))
+        self._assignment_mat = assignment_matrix
+        self._mitigation_mats = {}
+
+    @property
+    def settings(self) -> Dict:
+        """Return settings."""
+        return {"assignment_matrix": self._assignment_mat, "qubits": self._qubits}
+
+    def expectation_value(
+        self,
+        data: Counts,
+        diagonal: Union[Callable, dict, str, np.ndarray] = None,
+        qubits: Iterable[int] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> Tuple[float, float]:
+        r"""Compute the mitigated expectation value of a diagonal observable.
+
+        This computes the mitigated estimator of
+        :math:`\langle O \rangle = \mbox{Tr}[\rho. O]` of a diagonal observable
+        :math:`O = \sum_{x\in\{0, 1\}^n} O(x)|x\rangle\!\langle x|`.
+
+        Args:
+            data: Counts object
+            diagonal: Optional, the vector of diagonal values for summing the
+                      expectation value. If ``None`` the default value is
+                      :math:`[1, -1]^\otimes n`.
+            qubits: Optional, the measured physical qubits the count
+                    bitstrings correspond to. If None qubits are assumed to be
+                    :math:`[0, ..., n-1]`.
+            clbits: Optional, if not None marginalize counts to the specified bits.
+            shots: the number of shots.
+
+        Returns:
+            (float, float): the expectation value and an upper bound of the standard deviation.
+
+        Additional Information:
+            The diagonal observable :math:`O` is input using the ``diagonal`` kwarg as
+            a list or Numpy array :math:`[O(0), ..., O(2^n -1)]`. If no diagonal is specified
+            the diagonal of the Pauli operator
+            :math`O = \mbox{diag}(Z^{\otimes n}) = [1, -1]^{\otimes n}` is used.
+            The ``clbits`` kwarg is used to marginalize the input counts dictionary
+            over the specified bit-values, and the ``qubits`` kwarg is used to specify
+            which physical qubits these bit-values correspond to as
+            ``circuit.measure(qubits, clbits)``.
+        """
+
+        if qubits is None:
+            qubits = self._qubits
+        probs_vec, shots = counts_probability_vector(
+            data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits
+        )
+
+        # Get qubit mitigation matrix and mitigate probs
+        mit_mat = self.mitigation_matrix(qubits)
+
+        # Get operator coeffs
+        if diagonal is None:
+            diagonal = z_diagonal(2**self._num_qubits)
+        elif isinstance(diagonal, str):
+            diagonal = str2diag(diagonal)
+
+        # Apply transpose of mitigation matrix
+        coeffs = mit_mat.T.dot(diagonal)
+        expval = coeffs.dot(probs_vec)
+        stddev_upper_bound = self.stddev_upper_bound(shots)
+
+        return (expval, stddev_upper_bound)
+
+    def quasi_probabilities(
+        self,
+        data: Counts,
+        qubits: Optional[List[int]] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> QuasiDistribution:
+        """Compute mitigated quasi probabilities value.
+
+        Args:
+            data: counts object
+            qubits: qubits the count bitstrings correspond to.
+            clbits: Optional, marginalize counts to just these bits.
+            shots: Optional, the total number of shots, if None shots will
+                be calculated as the sum of all counts.
+
+        Returns:
+            QuasiDistribution: A dictionary containing pairs of [output, mean] where "output"
+                is the key in the dictionaries,
+                which is the length-N bitstring of a measured standard basis state,
+                and "mean" is the mean of non-zero quasi-probability estimates.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        probs_vec, calculated_shots = counts_probability_vector(
+            data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits
+        )
+        if shots is None:
+            shots = calculated_shots
+
+        # Get qubit mitigation matrix and mitigate probs
+        mit_mat = self.mitigation_matrix(qubits)
+
+        # Apply transpose of mitigation matrix
+        probs_vec = mit_mat.dot(probs_vec)
+        probs_dict = {}
+        for index, _ in enumerate(probs_vec):
+            probs_dict[index] = probs_vec[index]
+
+        quasi_dist = QuasiDistribution(
+            probs_dict, stddev_upper_bound=self.stddev_upper_bound(shots)
+        )
+
+        return quasi_dist
+
+    def mitigation_matrix(self, qubits: List[int] = None) -> np.ndarray:
+        r"""Return the readout mitigation matrix for the specified qubits.
+
+        The mitigation matrix :math:`A^{-1}` is defined as the inverse of the
+        :meth:`assignment_matrix` :math:`A`.
+
+        Args:
+            qubits: Optional, qubits being measured.
+
+        Returns:
+            np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        qubits = tuple(sorted(qubits))
+
+        # Check for cached mitigation matrix
+        # if not present compute
+        if qubits not in self._mitigation_mats:
+            marginal_matrix = self.assignment_matrix(qubits)
+            try:
+                mit_mat = np.linalg.inv(marginal_matrix)
+            except np.linalg.LinAlgError:
+                # Use pseudo-inverse if matrix is singular
+                mit_mat = np.linalg.pinv(marginal_matrix)
+            self._mitigation_mats[qubits] = mit_mat
+
+        return self._mitigation_mats[qubits]
+
+    def assignment_matrix(self, qubits: List[int] = None) -> np.ndarray:
+        r"""Return the readout assignment matrix for specified qubits.
+
+        The assignment matrix is the stochastic matrix :math:`A` which assigns
+        a noisy readout probability distribution to an ideal input
+        readout distribution: :math:`P(i|j) = \langle i | A | j \rangle`.
+
+        Args:
+            qubits: Optional, qubits being measured.
+
+        Returns:
+            np.ndarray: the assignment matrix A.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        if qubits == self._num_qubits:
+            return self._assignment_mat
+
+        if isinstance(qubits, int):
+            qubits = [qubits]
+
+        qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+        # Compute marginal matrix
+        axis = tuple(
+            self._num_qubits - 1 - i for i in set(range(self._num_qubits)).difference(qubit_indices)
+        )
+        num_qubits = len(qubits)
+
+        new_amat = np.zeros(2 * [2**num_qubits], dtype=float)
+        for i, col in enumerate(self._assignment_mat.T[self._keep_indexes(qubit_indices)]):
+            new_amat[i] = (
+                np.reshape(col, self._num_qubits * [2]).sum(axis=axis).reshape([2**num_qubits])
+            )
+        new_amat = new_amat.T
+        return new_amat
+
+    @staticmethod
+    def _keep_indexes(qubits):
+        indexes = [0]
+        for i in sorted(qubits):
+            indexes += [idx + (1 << i) for idx in indexes]
+        return indexes
+
+    def _compute_gamma(self):
+        """Compute gamma for N-qubit mitigation"""
+        mitmat = self.mitigation_matrix(qubits=self._qubits)
+        return np.max(np.sum(np.abs(mitmat), axis=0))
+
+    def stddev_upper_bound(self, shots: int):
+        """Return an upper bound on standard deviation of expval estimator.
+
+        Args:
+            shots: Number of shots used for expectation value measurement.
+
+        Returns:
+            float: the standard deviation upper bound.
+        """
+        gamma = self._compute_gamma()
+        return gamma / math.sqrt(shots)
+
+    @property
+    def qubits(self) -> Tuple[int]:
+        """The device qubits for this mitigator"""
+        return self._qubits
diff --git a/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py b/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py
new file mode 100644
index 0000000000..4e843dbbd4
--- /dev/null
+++ b/qiskit_experiments/data_processing/mitigation/local_readout_mitigator.py
@@ -0,0 +1,320 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Readout mitigator class based on the 1-qubit local tensored mitigation method
+"""
+
+
+import math
+from typing import Optional, List, Tuple, Iterable, Callable, Union, Dict
+import numpy as np
+
+from qiskit.exceptions import QiskitError
+from qiskit.result.distributions.quasi import QuasiDistribution
+from qiskit.result.counts import Counts
+from .base_readout_mitigator import BaseReadoutMitigator
+from .utils import counts_probability_vector, z_diagonal, str2diag
+
+
+class LocalReadoutMitigator(BaseReadoutMitigator):
+    """1-qubit tensor product readout error mitigator.
+
+    Mitigates :meth:`expectation_value` and :meth:`quasi_probabilities`.
+    The mitigator should either be calibrated using qiskit experiments,
+    or calculated directly from the backend properties.
+    This mitigation method should be used in case the readout errors of the qubits
+    are assumed to be uncorrelated. For *N* qubits there are *N* mitigation matrices,
+    each of size :math:`2 x 2` and the mitigation complexity is :math:`O(2^N)`,
+    so it is more efficient than the :class:`CorrelatedReadoutMitigator` class.
+    """
+
+    def __init__(
+        self,
+        assignment_matrices: Optional[List[np.ndarray]] = None,
+        qubits: Optional[Iterable[int]] = None,
+        backend=None,
+    ):
+        """Initialize a LocalReadoutMitigator
+
+        Args:
+            assignment_matrices: Optional, list of single-qubit readout error assignment matrices.
+            qubits: Optional, the measured physical qubits for mitigation.
+            backend: Optional, backend name.
+
+        Raises:
+            QiskitError: matrices sizes do not agree with number of qubits
+        """
+        if assignment_matrices is None:
+            assignment_matrices = self._from_backend(backend, qubits)
+        else:
+            assignment_matrices = [np.asarray(amat, dtype=float) for amat in assignment_matrices]
+        for amat in assignment_matrices:
+            if np.any(amat < 0) or not np.allclose(np.sum(amat, axis=0), 1):
+                raise QiskitError(
+                    "Assignment matrix columns must be valid probability distributions"
+                )
+        if qubits is None:
+            self._num_qubits = len(assignment_matrices)
+            self._qubits = range(self._num_qubits)
+        else:
+            if len(qubits) != len(assignment_matrices):
+                raise QiskitError(
+                    f"The number of given qubits ({len(qubits)}) is different than the number of qubits "
+                    f"inferred from the matrices ({len(assignment_matrices)})"
+                )
+            self._qubits = qubits
+            self._num_qubits = len(self._qubits)
+
+        self._qubit_index = dict(zip(self._qubits, range(self._num_qubits)))
+        self._assignment_mats = assignment_matrices
+        self._mitigation_mats = np.zeros([self._num_qubits, 2, 2], dtype=float)
+        self._gammas = np.zeros(self._num_qubits, dtype=float)
+
+        for i in range(self._num_qubits):
+            mat = self._assignment_mats[i]
+            # Compute Gamma values
+            error0 = mat[1, 0]
+            error1 = mat[0, 1]
+            self._gammas[i] = (1 + abs(error0 - error1)) / (1 - error0 - error1)
+            # Compute inverse mitigation matrix
+            try:
+                ainv = np.linalg.inv(mat)
+            except np.linalg.LinAlgError:
+                ainv = np.linalg.pinv(mat)
+            self._mitigation_mats[i] = ainv
+
+    @property
+    def settings(self) -> Dict:
+        """Return settings."""
+        return {"assignment_matrices": self._assignment_mats, "qubits": self._qubits}
+
+    def expectation_value(
+        self,
+        data: Counts,
+        diagonal: Union[Callable, dict, str, np.ndarray] = None,
+        qubits: Iterable[int] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> Tuple[float, float]:
+        r"""Compute the mitigated expectation value of a diagonal observable.
+
+        This computes the mitigated estimator of
+        :math:`\langle O \rangle = \mbox{Tr}[\rho. O]` of a diagonal observable
+        :math:`O = \sum_{x\in\{0, 1\}^n} O(x)|x\rangle\!\langle x|`.
+
+        Args:
+            data: Counts object
+            diagonal: Optional, the vector of diagonal values for summing the
+                      expectation value. If ``None`` the default value is
+                      :math:`[1, -1]^\otimes n`.
+            qubits: Optional, the measured physical qubits the count
+                    bitstrings correspond to. If None qubits are assumed to be
+                    :math:`[0, ..., n-1]`.
+            clbits: Optional, if not None marginalize counts to the specified bits.
+            shots: the number of shots.
+
+        Returns:
+            (float, float): the expectation value and an upper bound of the standard deviation.
+
+        Additional Information:
+            The diagonal observable :math:`O` is input using the ``diagonal`` kwarg as
+            a list or Numpy array :math:`[O(0), ..., O(2^n -1)]`. If no diagonal is specified
+            the diagonal of the Pauli operator
+            :math`O = \mbox{diag}(Z^{\otimes n}) = [1, -1]^{\otimes n}` is used.
+            The ``clbits`` kwarg is used to marginalize the input counts dictionary
+            over the specified bit-values, and the ``qubits`` kwarg is used to specify
+            which physical qubits these bit-values correspond to as
+            ``circuit.measure(qubits, clbits)``.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        num_qubits = len(qubits)
+        probs_vec, shots = counts_probability_vector(
+            data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits
+        )
+
+        # Get qubit mitigation matrix and mitigate probs
+        qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+        ainvs = self._mitigation_mats[qubit_indices]
+
+        # Get operator coeffs
+        if diagonal is None:
+            diagonal = z_diagonal(2**num_qubits)
+        elif isinstance(diagonal, str):
+            diagonal = str2diag(diagonal)
+
+        # Apply transpose of mitigation matrix
+        coeffs = np.reshape(diagonal, num_qubits * [2])
+        einsum_args = [coeffs, list(range(num_qubits))]
+        for i, ainv in enumerate(reversed(ainvs)):
+            einsum_args += [ainv.T, [num_qubits + i, i]]
+        einsum_args += [list(range(num_qubits, 2 * num_qubits))]
+        coeffs = np.einsum(*einsum_args).ravel()
+
+        expval = coeffs.dot(probs_vec)
+        stddev_upper_bound = self.stddev_upper_bound(shots, qubits)
+
+        return (expval, stddev_upper_bound)
+
+    def quasi_probabilities(
+        self,
+        data: Counts,
+        qubits: Optional[List[int]] = None,
+        clbits: Optional[List[int]] = None,
+        shots: Optional[int] = None,
+    ) -> QuasiDistribution:
+        """Compute mitigated quasi probabilities value.
+
+        Args:
+            data: counts object
+            qubits: qubits the count bitstrings correspond to.
+            clbits: Optional, marginalize counts to just these bits.
+            shots: Optional, the total number of shots, if None shots will
+                be calculated as the sum of all counts.
+
+        Returns:
+            QuasiDistribution: A dictionary containing pairs of [output, mean] where "output"
+                is the key in the dictionaries,
+                which is the length-N bitstring of a measured standard basis state,
+                and "mean" is the mean of non-zero quasi-probability estimates.
+
+        Raises:
+            QiskitError: if qubit and clbit kwargs are not valid.
+        """
+        if qubits is None:
+            qubits = self._qubits
+
+        num_qubits = len(qubits)
+
+        probs_vec, calculated_shots = counts_probability_vector(
+            data, qubit_index=self._qubit_index, clbits=clbits, qubits=qubits
+        )
+        if shots is None:
+            shots = calculated_shots
+
+        # Get qubit mitigation matrix and mitigate probs
+        qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+        ainvs = self._mitigation_mats[qubit_indices]
+
+        # Apply transpose of mitigation matrix
+        prob_tens = np.reshape(probs_vec, num_qubits * [2])
+        einsum_args = [prob_tens, list(range(num_qubits))]
+        for i, ainv in enumerate(reversed(ainvs)):
+            einsum_args += [ainv, [num_qubits + i, i]]
+        einsum_args += [list(range(num_qubits, 2 * num_qubits))]
+        probs_vec = np.einsum(*einsum_args).ravel()
+
+        probs_dict = {}
+        for index, _ in enumerate(probs_vec):
+            probs_dict[index] = probs_vec[index]
+
+        quasi_dist = QuasiDistribution(
+            probs_dict, shots=shots, stddev_upper_bound=self.stddev_upper_bound(shots, qubits)
+        )
+        return quasi_dist
+
+    def mitigation_matrix(self, qubits: Optional[Union[List[int], int]] = None) -> np.ndarray:
+        r"""Return the measurement mitigation matrix for the specified qubits.
+
+        The mitigation matrix :math:`A^{-1}` is defined as the inverse of the
+        :meth:`assignment_matrix` :math:`A`.
+
+        Args:
+            qubits: Optional, qubits being measured for operator expval.
+                    if a single int is given, it is assumed to be the index
+                    of the qubit in self._qubits
+
+        Returns:
+            np.ndarray: the measurement error mitigation matrix :math:`A^{-1}`.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        if isinstance(qubits, int):
+            qubits = [self._qubits[qubits]]
+        qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+        mat = self._mitigation_mats[qubit_indices[0]]
+        for i in qubit_indices[1:]:
+            mat = np.kron(self._mitigation_mats[i], mat)
+        return mat
+
+    def assignment_matrix(self, qubits: List[int] = None) -> np.ndarray:
+        r"""Return the measurement assignment matrix for specified qubits.
+
+        The assignment matrix is the stochastic matrix :math:`A` which assigns
+        a noisy measurement probability distribution to an ideal input
+        measurement distribution: :math:`P(i|j) = \langle i | A | j \rangle`.
+
+        Args:
+            qubits: Optional, qubits being measured for operator expval.
+
+        Returns:
+            np.ndarray: the assignment matrix A.
+        """
+        if qubits is None:
+            qubits = self._qubits
+        if isinstance(qubits, int):
+            qubits = [qubits]
+        qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+        mat = self._assignment_mats[qubit_indices[0]]
+        for i in qubit_indices[1:]:
+            mat = np.kron(self._assignment_mats[i], mat)
+        return mat
+
+    def _compute_gamma(self, qubits=None):
+        """Compute gamma for N-qubit mitigation"""
+        if qubits is None:
+            gammas = self._gammas
+        else:
+            qubit_indices = [self._qubit_index[qubit] for qubit in qubits]
+            gammas = self._gammas[qubit_indices]
+        return np.prod(gammas)
+
+    def stddev_upper_bound(self, shots: int, qubits: List[int] = None):
+        """Return an upper bound on standard deviation of expval estimator.
+
+        Args:
+            shots: Number of shots used for expectation value measurement.
+            qubits: qubits being measured for operator expval.
+
+        Returns:
+            float: the standard deviation upper bound.
+        """
+        gamma = self._compute_gamma(qubits=qubits)
+        return gamma / math.sqrt(shots)
+
+    def _from_backend(self, backend, qubits):
+        """Calculates amats from backend properties readout_error"""
+        backend_qubits = backend.properties().qubits
+        if qubits is not None:
+            if any(qubit >= len(backend_qubits) for qubit in qubits):
+                raise QiskitError("The chosen backend does not contain the specified qubits.")
+            reduced_backend_qubits = [backend_qubits[i] for i in qubits]
+            backend_qubits = reduced_backend_qubits
+        num_qubits = len(backend_qubits)
+
+        amats = np.zeros([num_qubits, 2, 2], dtype=float)
+
+        for qubit_idx, qubit_prop in enumerate(backend_qubits):
+            for prop in qubit_prop:
+                if prop.name == "prob_meas0_prep1":
+                    (amats[qubit_idx])[0, 1] = prop.value
+                    (amats[qubit_idx])[1, 1] = 1 - prop.value
+                if prop.name == "prob_meas1_prep0":
+                    (amats[qubit_idx])[1, 0] = prop.value
+                    (amats[qubit_idx])[0, 0] = 1 - prop.value
+
+        return amats
+
+    @property
+    def qubits(self) -> Tuple[int]:
+        """The device qubits for this mitigator"""
+        return self._qubits
diff --git a/qiskit_experiments/data_processing/mitigation/utils.py b/qiskit_experiments/data_processing/mitigation/utils.py
new file mode 100644
index 0000000000..6f4b9d7871
--- /dev/null
+++ b/qiskit_experiments/data_processing/mitigation/utils.py
@@ -0,0 +1,163 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2021
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+"""
+Readout mitigation data handling utils
+
+.. note::
+
+    None of these functions are intended for public use.
+"""
+
+import logging
+import math
+from typing import Optional, List, Tuple, Dict
+import numpy as np
+
+from qiskit.exceptions import QiskitError
+from qiskit.result import Counts, marginal_counts
+
+logger = logging.getLogger(__name__)
+
+
+def z_diagonal(dim, dtype=float):
+    r"""Return the diagonal for the operator :math:`Z^\otimes n`"""
+    parity = np.zeros(dim, dtype=dtype)
+    for i in range(dim):
+        parity[i] = bin(i)[2:].count("1")
+    return (-1) ** np.mod(parity, 2)
+
+
+def expval_with_stddev(coeffs: np.ndarray, probs: np.ndarray, shots: int) -> Tuple[float, float]:
+    """Compute expectation value and standard deviation.
+    Args:
+        coeffs: array of diagonal operator coefficients.
+        probs: array of measurement probabilities.
+        shots: total number of shots to obtain probabilities.
+    Returns:
+        tuple: (expval, stddev) expectation value and standard deviation.
+    """
+    # Compute expval
+    expval = coeffs.dot(probs)
+
+    # Compute variance
+    sq_expval = (coeffs**2).dot(probs)
+    variance = (sq_expval - expval**2) / shots
+
+    # Compute standard deviation
+    if variance < 0 and not np.isclose(variance, 0):
+        logger.warning(
+            "Encountered a negative variance in expectation value calculation."
+            "(%f). Setting standard deviation of result to 0.",
+            variance,
+        )
+    calc_stddev = math.sqrt(variance) if variance > 0 else 0.0
+    return [expval, calc_stddev]
+
+
+def stddev(probs, shots):
+    """Calculate stddev dict"""
+    ret = {}
+    for key, prob in probs.items():
+        std_err = math.sqrt(prob * (1 - prob) / shots)
+        ret[key] = std_err
+    return ret
+
+
+def str2diag(string):
+    """Transform diagonal from a string to a numpy array"""
+    chars = {
+        "I": np.array([1, 1], dtype=float),
+        "Z": np.array([1, -1], dtype=float),
+        "0": np.array([1, 0], dtype=float),
+        "1": np.array([0, 1], dtype=float),
+    }
+    ret = np.array([1], dtype=float)
+    for i in reversed(string):
+        if i not in chars:
+            raise QiskitError(f"Invalid diagonal string character {i}")
+        ret = np.kron(chars[i], ret)
+    return ret
+
+
+def counts_to_vector(counts: Counts, num_qubits: int) -> Tuple[np.ndarray, int]:
+    """Transforms Counts to a probability vector"""
+    vec = np.zeros(2**num_qubits, dtype=float)
+    shots = 0
+    for key, val in counts.items():
+        shots += val
+        vec[int(key, 2)] = val
+    vec /= shots
+    return vec, shots
+
+
+def remap_qubits(
+    vec: np.ndarray, num_qubits: int, qubits: Optional[List[int]] = None
+) -> np.ndarray:
+    """Remapping the qubits"""
+    if qubits is not None:
+        if len(qubits) != num_qubits:
+            raise QiskitError("Num qubits does not match vector length.")
+        axes = [num_qubits - 1 - i for i in reversed(np.argsort(qubits))]
+        vec = np.reshape(vec, num_qubits * [2]).transpose(axes).reshape(vec.shape)
+    return vec
+
+
+def marginalize_counts(
+    counts: Counts,
+    qubit_index: Dict[int, int],
+    qubits: Optional[List[int]] = None,
+    clbits: Optional[List[int]] = None,
+) -> np.ndarray:
+    """Marginalization of the Counts. Verify that number of clbits equals to the number of qubits."""
+    if clbits is not None:
+        qubits_len = len(qubits) if not qubits is None else 0
+        clbits_len = len(clbits) if not clbits is None else 0
+        if clbits_len not in (0, qubits_len):
+            raise QiskitError(
+                f"Num qubits ({qubits_len}) does not match number of clbits ({clbits_len})."
+            )
+        counts = marginal_counts(counts, clbits)
+    if clbits is None and qubits is not None:
+        clbits = [qubit_index[qubit] for qubit in qubits]
+        counts = marginal_counts(counts, clbits)
+    return counts
+
+
+def counts_probability_vector(
+    counts: Counts,
+    qubit_index: Dict[int, int],
+    qubits: Optional[List[int]] = None,
+    clbits: Optional[List[int]] = None,
+) -> Tuple[np.ndarray, int]:
+    """Compute a probability vector for all count outcomes.
+
+    Args:
+        counts: counts object
+        qubit_index: For each qubit, its index in the mitigator qubits list
+        qubits: qubits the count bitstrings correspond to.
+        clbits: Optional, marginalize counts to just these bits.
+
+    Raises:
+        QiskitError: if qubits and clbits kwargs are not valid.
+
+    Returns:
+        np.ndarray: a probability vector for all count outcomes.
+        int: Number of shots in the counts
+    """
+    counts = marginalize_counts(counts, qubit_index, qubits, clbits)
+    if qubits is not None:
+        num_qubits = len(qubits)
+    else:
+        num_qubits = len(qubit_index.keys())
+    vec, shots = counts_to_vector(counts, num_qubits)
+    vec = remap_qubits(vec, num_qubits, qubits)
+    return vec, shots
diff --git a/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py b/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py
index 73dece5171..796c375d9f 100644
--- a/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py
+++ b/qiskit_experiments/library/characterization/analysis/correlated_readout_error_analysis.py
@@ -15,7 +15,7 @@
 from typing import List, Tuple
 import numpy as np
 import matplotlib.pyplot as plt
-from qiskit.result import CorrelatedReadoutMitigator
+from qiskit_experiments.data_processing import CorrelatedReadoutMitigator
 from qiskit_experiments.framework import ExperimentData
 from qiskit_experiments.framework.matplotlib import get_non_gui_ax
 from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options
diff --git a/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py b/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py
index 0bb3263496..ee8917fb0e 100644
--- a/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py
+++ b/qiskit_experiments/library/characterization/analysis/local_readout_error_analysis.py
@@ -15,8 +15,8 @@
 from typing import List, Tuple
 import numpy as np
 import matplotlib.pyplot as plt
-from qiskit.result import LocalReadoutMitigator
 from qiskit.result import marginal_distribution
+from qiskit_experiments.data_processing import LocalReadoutMitigator
 from qiskit_experiments.framework import ExperimentData
 from qiskit_experiments.framework.matplotlib import get_non_gui_ax
 from qiskit_experiments.framework import BaseAnalysis, AnalysisResultData, Options
diff --git a/qiskit_experiments/library/tomography/basis/pauli_basis.py b/qiskit_experiments/library/tomography/basis/pauli_basis.py
index 1d908ef322..082da12333 100644
--- a/qiskit_experiments/library/tomography/basis/pauli_basis.py
+++ b/qiskit_experiments/library/tomography/basis/pauli_basis.py
@@ -16,9 +16,9 @@
 import numpy as np
 from qiskit.circuit import QuantumCircuit
 from qiskit.circuit.library import HGate, XGate, ZGate, SGate, SdgGate
-from qiskit.result import LocalReadoutMitigator
 from qiskit.quantum_info import DensityMatrix
 from qiskit.exceptions import QiskitError
+from qiskit_experiments.data_processing import LocalReadoutMitigator
 from .local_basis import LocalMeasurementBasis, LocalPreparationBasis
 
 
diff --git a/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml b/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml
new file mode 100644
index 0000000000..61a06234a4
--- /dev/null
+++ b/releasenotes/notes/adopt-mitigators-d7ad0e2f3cd2fa57.yaml
@@ -0,0 +1,18 @@
+---
+features:
+  - |
+    New :class:`~.LocalReadoutMitigator` and
+    :class:`~.CorrelatedReadoutMitigator` classes have been added. These
+    classes were moved directly from Qiskit which deprecated them in Qiskit
+    1.3. They provide utility methods for applying readout error mitigation and
+    integrate with the readout error mitigation experiments
+    :class:`~.LocalReadoutError` and :class:`~.CorrelatedReadoutError`.
+upgrade:
+  - |
+    The readout error mitigation experiments :class:`~.LocalReadoutError` and
+    :class:`~.CorrelatedReadoutError` have been updated to generate instances
+    of the new :class:`~.LocalReadoutMitigator` and
+    :class:`~.CorrelatedReadoutMitigator` classes. The experiments should
+    continue to work as before, but any code that was using, for example,
+    ``isinstance()`` to check object type would need to be updated to check
+    against the Qiskit Experiments classes instead of the old Qiskit classes.
diff --git a/test/data_processing/test_mitigators.py b/test/data_processing/test_mitigators.py
new file mode 100644
index 0000000000..b75f7fde12
--- /dev/null
+++ b/test/data_processing/test_mitigators.py
@@ -0,0 +1,477 @@
+# This code is part of Qiskit.
+#
+# (C) Copyright IBM 2017, 2024.
+#
+# This code is licensed under the Apache License, Version 2.0. You may
+# obtain a copy of this license in the LICENSE.txt file in the root directory
+# of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
+#
+# Any modifications or derivative works of this code must retain this
+# copyright notice, and modified files need to carry a notice indicating
+# that they have been altered from the originals.
+# pylint: disable=invalid-name
+
+"""Tests for error mitigation routines."""
+
+import unittest
+from collections import Counter
+from test.base import QiskitExperimentsTestCase
+
+import numpy as np
+
+from qiskit import QiskitError, QuantumCircuit
+from qiskit.quantum_info import Statevector
+from qiskit.quantum_info.operators.predicates import matrix_equal
+from qiskit.result import Counts
+from qiskit.result.utils import marginal_counts
+from qiskit.providers.fake_provider import Fake5QV1
+
+from qiskit_experiments.data_processing import (
+    CorrelatedReadoutMitigator,
+    LocalReadoutMitigator,
+)
+from qiskit_experiments.data_processing.mitigation import (
+    counts_probability_vector,
+    expval_with_stddev,
+    stddev,
+    str2diag,
+)
+
+
+class TestReadoutMitigation(QiskitExperimentsTestCase):
+    """Tests for correlated and local readout mitigation."""
+
+    @classmethod
+    def setUpClass(cls):
+        super().setUpClass()
+        cls.rng = np.random.default_rng(42)
+
+    @staticmethod
+    def compare_results(res1, res2):
+        """Compare the results between two runs"""
+        res1_total_shots = sum(res1.values())
+        res2_total_shots = sum(res2.values())
+        keys = set(res1.keys()).union(set(res2.keys()))
+        total = 0
+        for key in keys:
+            val1 = res1.get(key, 0) / res1_total_shots
+            val2 = res2.get(key, 0) / res2_total_shots
+            total += abs(val1 - val2) ** 2
+        return total
+
+    @staticmethod
+    def mitigators(assignment_matrices, qubits=None):
+        """Generates the mitigators to test for given assignment matrices"""
+        full_assignment_matrix = assignment_matrices[0]
+        for m in assignment_matrices[1:]:
+            full_assignment_matrix = np.kron(full_assignment_matrix, m)
+        CRM = CorrelatedReadoutMitigator(full_assignment_matrix, qubits)
+        LRM = LocalReadoutMitigator(assignment_matrices, qubits)
+        mitigators = [CRM, LRM]
+        return mitigators
+
+    @staticmethod
+    def simulate_circuit(circuit, assignment_matrix, num_qubits, shots=1024):
+        """Simulates the given circuit under the given readout noise"""
+        probs = Statevector.from_instruction(circuit).probabilities()
+        noisy_probs = assignment_matrix @ probs
+        labels = [bin(a)[2:].zfill(num_qubits) for a in range(2**num_qubits)]
+        results = TestReadoutMitigation.rng.choice(labels, size=shots, p=noisy_probs)
+        return Counts(dict(Counter(results)))
+
+    @staticmethod
+    def ghz_3_circuit():
+        """A 3-qubit circuit generating |000>+|111>"""
+        c = QuantumCircuit(3)
+        c.h(0)
+        c.cx(0, 1)
+        c.cx(1, 2)
+        return (c, "ghz_3_ciruit", 3)
+
+    @staticmethod
+    def first_qubit_h_3_circuit():
+        """A 3-qubit circuit generating |000>+|001>"""
+        c = QuantumCircuit(3)
+        c.h(0)
+        return (c, "first_qubit_h_3_circuit", 3)
+
+    @staticmethod
+    def assignment_matrices():
+        """A 3-qubit readout noise assignment matrices"""
+        return LocalReadoutMitigator(backend=Fake5QV1())._assignment_mats[0:3]
+
+    @staticmethod
+    def counts_data(circuit, assignment_matrices, shots=1024):
+        """Generates count data for the noisy and noiseless versions of the circuit simulation"""
+        full_assignment_matrix = assignment_matrices[0]
+        for m in assignment_matrices[1:]:
+            full_assignment_matrix = np.kron(full_assignment_matrix, m)
+        num_qubits = len(assignment_matrices)
+        ideal_assignment_matrix = np.eye(2**num_qubits)
+        counts_ideal = TestReadoutMitigation.simulate_circuit(
+            circuit, ideal_assignment_matrix, num_qubits, shots
+        )
+        counts_noise = TestReadoutMitigation.simulate_circuit(
+            circuit, full_assignment_matrix, num_qubits, shots
+        )
+        probs_noise = {key: value / shots for key, value in counts_noise.items()}
+        return counts_ideal, counts_noise, probs_noise
+
+    def test_mitigation_improvement(self):
+        """Test whether readout mitigation led to more accurate results"""
+        shots = 1024
+        with self.assertWarns(DeprecationWarning):
+            # TODO self.assignment_matrices calls LocalReadoutMitigator,
+            #  which only supports BackendV1 at the moment:
+            #  https://github.com/Qiskit/qiskit/issues/12832
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices)
+        circuit, circuit_name, num_qubits = self.ghz_3_circuit()
+        counts_ideal, counts_noise, probs_noise = self.counts_data(
+            circuit, assignment_matrices, shots
+        )
+        unmitigated_error = self.compare_results(counts_ideal, counts_noise)
+        unmitigated_stddev = stddev(probs_noise, shots)
+
+        for mitigator in mitigators:
+            mitigated_quasi_probs = mitigator.quasi_probabilities(counts_noise)
+            mitigated_probs = (
+                mitigated_quasi_probs.nearest_probability_distribution().binary_probabilities(
+                    num_bits=num_qubits
+                )
+            )
+            mitigated_error = self.compare_results(counts_ideal, mitigated_probs)
+            self.assertLess(
+                mitigated_error,
+                unmitigated_error * 0.8,
+                f"Mitigator {mitigator} did not improve circuit {circuit_name} measurements",
+            )
+            mitigated_stddev_upper_bound = mitigated_quasi_probs._stddev_upper_bound
+            max_unmitigated_stddev = max(unmitigated_stddev.values())
+            self.assertGreaterEqual(
+                mitigated_stddev_upper_bound,
+                max_unmitigated_stddev,
+                f"Mitigator {mitigator} on circuit {circuit_name} gave stddev upper bound "
+                f"{mitigated_stddev_upper_bound} while unmitigated stddev maximum is "
+                f"{max_unmitigated_stddev}",
+            )
+
+    def test_expectation_improvement(self):
+        """Test whether readout mitigation led to more accurate results
+        and that its standard deviation is increased"""
+        shots = 1024
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices)
+        num_qubits = len(assignment_matrices)
+        diagonals = []
+        diagonals.append("IZ0")
+        diagonals.append("101")
+        diagonals.append("IZZ")
+        qubit_index = {i: i for i in range(num_qubits)}
+        circuit, circuit_name, num_qubits = self.ghz_3_circuit()
+        counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots)
+        probs_ideal, _ = counts_probability_vector(counts_ideal, qubit_index=qubit_index)
+        probs_noise, _ = counts_probability_vector(counts_noise, qubit_index=qubit_index)
+        for diagonal in diagonals:
+            if isinstance(diagonal, str):
+                diagonal = str2diag(diagonal)
+            unmitigated_expectation, unmitigated_stddev = expval_with_stddev(
+                diagonal, probs_noise, shots=counts_noise.shots()
+            )
+            ideal_expectation = np.dot(probs_ideal, diagonal)
+            unmitigated_error = np.abs(ideal_expectation - unmitigated_expectation)
+            for mitigator in mitigators:
+                mitigated_expectation, mitigated_stddev = mitigator.expectation_value(
+                    counts_noise, diagonal
+                )
+                mitigated_error = np.abs(ideal_expectation - mitigated_expectation)
+                self.assertLess(
+                    mitigated_error,
+                    unmitigated_error,
+                    f"Mitigator {mitigator} did not improve circuit {circuit_name} expectation "
+                    f"computation for diagonal {diagonal} ideal: {ideal_expectation}, unmitigated:"
+                    f" {unmitigated_expectation} mitigated: {mitigated_expectation}",
+                )
+                self.assertGreaterEqual(
+                    mitigated_stddev,
+                    unmitigated_stddev,
+                    f"Mitigator {mitigator} did not increase circuit {circuit_name} the"
+                    f" standard deviation",
+                )
+
+    def test_clbits_parameter(self):
+        """Test whether the clbits parameter is handled correctly"""
+        shots = 10000
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices)
+        circuit, _, _ = self.first_qubit_h_3_circuit()
+        counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots)
+        counts_ideal_12 = marginal_counts(counts_ideal, [1, 2])
+        counts_ideal_02 = marginal_counts(counts_ideal, [0, 2])
+
+        for mitigator in mitigators:
+            mitigated_probs_12 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[1, 2], clbits=[1, 2])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=2)
+            )
+            mitigated_error = self.compare_results(counts_ideal_12, mitigated_probs_12)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly marginalize for qubits 1,2",
+            )
+
+            mitigated_probs_02 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[0, 2], clbits=[0, 2])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=2)
+            )
+            mitigated_error = self.compare_results(counts_ideal_02, mitigated_probs_02)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly marginalize for qubits 0,2",
+            )
+
+    def test_qubits_parameter(self):
+        """Test whether the qubits parameter is handled correctly"""
+        shots = 10000
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices)
+        circuit, _, _ = self.first_qubit_h_3_circuit()
+        counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots)
+        counts_ideal_012 = counts_ideal
+        counts_ideal_210 = Counts({"000": counts_ideal["000"], "100": counts_ideal["001"]})
+        counts_ideal_102 = Counts({"000": counts_ideal["000"], "010": counts_ideal["001"]})
+
+        for mitigator in mitigators:
+            mitigated_probs_012 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[0, 1, 2])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=3)
+            )
+            mitigated_error = self.compare_results(counts_ideal_012, mitigated_probs_012)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly handle qubit order 0, 1, 2",
+            )
+
+            mitigated_probs_210 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[2, 1, 0])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=3)
+            )
+            mitigated_error = self.compare_results(counts_ideal_210, mitigated_probs_210)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly handle qubit order 2, 1, 0",
+            )
+
+            mitigated_probs_102 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[1, 0, 2])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=3)
+            )
+            mitigated_error = self.compare_results(counts_ideal_102, mitigated_probs_102)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                "Mitigator {mitigator} did not correctly handle qubit order 1, 0, 2",
+            )
+
+    def test_repeated_qubits_parameter(self):
+        """Tests the order of mitigated qubits."""
+        shots = 10000
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices, qubits=[0, 1, 2])
+        circuit, _, _ = self.first_qubit_h_3_circuit()
+        counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots)
+        counts_ideal_012 = counts_ideal
+        counts_ideal_210 = Counts({"000": counts_ideal["000"], "100": counts_ideal["001"]})
+
+        for mitigator in mitigators:
+            mitigated_probs_210 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[2, 1, 0])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=3)
+            )
+            mitigated_error = self.compare_results(counts_ideal_210, mitigated_probs_210)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly handle qubit order 2,1,0",
+            )
+
+            # checking qubit order 2,1,0 should not "overwrite" the default 0,1,2
+            mitigated_probs_012 = (
+                mitigator.quasi_probabilities(counts_noise)
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=3)
+            )
+            mitigated_error = self.compare_results(counts_ideal_012, mitigated_probs_012)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly handle qubit order 0,1,2 "
+                f"(the expected default)",
+            )
+
+    def test_qubits_subset_parameter(self):
+        """Tests mitigation on a subset of the initial set of qubits."""
+
+        shots = 10000
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices, qubits=[2, 4, 6])
+        circuit, _, _ = self.first_qubit_h_3_circuit()
+        counts_ideal, counts_noise, _ = self.counts_data(circuit, assignment_matrices, shots)
+        counts_ideal_2 = marginal_counts(counts_ideal, [0])
+        counts_ideal_6 = marginal_counts(counts_ideal, [2])
+
+        for mitigator in mitigators:
+            mitigated_probs_2 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[2])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=1)
+            )
+            mitigated_error = self.compare_results(counts_ideal_2, mitigated_probs_2)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                "Mitigator {mitigator} did not correctly handle qubit subset",
+            )
+
+            mitigated_probs_6 = (
+                mitigator.quasi_probabilities(counts_noise, qubits=[6])
+                .nearest_probability_distribution()
+                .binary_probabilities(num_bits=1)
+            )
+            mitigated_error = self.compare_results(counts_ideal_6, mitigated_probs_6)
+            self.assertLess(
+                mitigated_error,
+                0.001,
+                f"Mitigator {mitigator} did not correctly handle qubit subset",
+            )
+            diagonal = str2diag("ZZ")
+            ideal_expectation = 0
+            mitigated_expectation, _ = mitigator.expectation_value(
+                counts_noise, diagonal, qubits=[2, 6]
+            )
+            mitigated_error = np.abs(ideal_expectation - mitigated_expectation)
+            self.assertLess(
+                mitigated_error,
+                0.1,
+                f"Mitigator {mitigator} did not improve circuit expectation",
+            )
+
+    def test_from_backend(self):
+        """Test whether a local mitigator can be created directly from backend properties"""
+        with self.assertWarns(DeprecationWarning):
+            backend = Fake5QV1()
+        num_qubits = len(backend.properties().qubits)
+        probs = TestReadoutMitigation.rng.random((num_qubits, 2))
+        for qubit_idx, qubit_prop in enumerate(backend.properties().qubits):
+            for prop in qubit_prop:
+                if prop.name == "prob_meas1_prep0":
+                    prop.value = probs[qubit_idx][0]
+                if prop.name == "prob_meas0_prep1":
+                    prop.value = probs[qubit_idx][1]
+        LRM_from_backend = LocalReadoutMitigator(backend=backend)
+
+        mats = []
+        for qubit_idx in range(num_qubits):
+            mat = np.array(
+                [
+                    [1 - probs[qubit_idx][0], probs[qubit_idx][1]],
+                    [probs[qubit_idx][0], 1 - probs[qubit_idx][1]],
+                ]
+            )
+            mats.append(mat)
+        LRM_from_matrices = LocalReadoutMitigator(assignment_matrices=mats)
+        self.assertTrue(
+            matrix_equal(
+                LRM_from_backend.assignment_matrix(), LRM_from_matrices.assignment_matrix()
+            )
+        )
+
+    def test_error_handling(self):
+        """Test that the assignment matrices are valid."""
+        bad_matrix_A = np.array([[-0.3, 1], [1.3, 0]])  # negative indices
+        bad_matrix_B = np.array([[0.2, 1], [0.7, 0]])  # columns not summing to 1
+        good_matrix_A = np.array([[0.2, 1], [0.8, 0]])
+        for bad_matrix in [bad_matrix_A, bad_matrix_B]:
+            with self.assertRaises(QiskitError) as cm:
+                CorrelatedReadoutMitigator(bad_matrix)
+            self.assertEqual(
+                cm.exception.message,
+                "Assignment matrix columns must be valid probability distributions",
+            )
+
+        with self.assertRaises(QiskitError) as cm:
+            amats = [good_matrix_A, bad_matrix_A]
+            LocalReadoutMitigator(amats)
+        self.assertEqual(
+            cm.exception.message,
+            "Assignment matrix columns must be valid probability distributions",
+        )
+
+        with self.assertRaises(QiskitError) as cm:
+            amats = [bad_matrix_B, good_matrix_A]
+            LocalReadoutMitigator(amats)
+        self.assertEqual(
+            cm.exception.message,
+            "Assignment matrix columns must be valid probability distributions",
+        )
+
+    def test_expectation_value_endian(self):
+        """Test that endian for expval is little."""
+        with self.assertWarns(DeprecationWarning):
+            assignment_matrices = self.assignment_matrices()
+        mitigators = self.mitigators(assignment_matrices)
+        counts = Counts({"10": 3, "11": 24, "00": 74, "01": 923})
+        for mitigator in mitigators:
+            expval, _ = mitigator.expectation_value(counts, diagonal="IZ", qubits=[0, 1])
+            self.assertAlmostEqual(expval, -1.0, places=0)
+
+    def test_quasi_probabilities_shots_passing(self):
+        """Test output of LocalReadoutMitigator.quasi_probabilities
+
+        We require the number of shots to be set in the output.
+        """
+        mitigator = LocalReadoutMitigator([np.array([[0.9, 0.1], [0.1, 0.9]])], qubits=[0])
+        counts = Counts({"10": 3, "11": 24, "00": 74, "01": 923})
+        quasi_dist = mitigator.quasi_probabilities(counts)
+        self.assertEqual(quasi_dist.shots, sum(counts.values()))
+
+        # custom number of shots
+        quasi_dist = mitigator.quasi_probabilities(counts, shots=1025)
+        self.assertEqual(quasi_dist.shots, 1025)
+
+
+class TestLocalReadoutMitigation(QiskitExperimentsTestCase):
+    """Tests specific to the local readout mitigator"""
+
+    def test_assignment_matrix(self):
+        """Tests that the local mitigator generates the full assignment matrix correctly"""
+        qubits = [7, 2, 3]
+        with self.assertWarns(DeprecationWarning):
+            backend = Fake5QV1()
+        assignment_matrices = LocalReadoutMitigator(backend=backend)._assignment_mats[0:3]
+        expected_assignment_matrix = np.kron(
+            np.kron(assignment_matrices[2], assignment_matrices[1]), assignment_matrices[0]
+        )
+        expected_mitigation_matrix = np.linalg.inv(expected_assignment_matrix)
+        LRM = LocalReadoutMitigator(assignment_matrices, qubits)
+        self.assertTrue(matrix_equal(expected_mitigation_matrix, LRM.mitigation_matrix()))
+        self.assertTrue(matrix_equal(expected_assignment_matrix, LRM.assignment_matrix()))
+
+
+if __name__ == "__main__":
+    unittest.main()
diff --git a/test/library/tomography/test_process_tomography.py b/test/library/tomography/test_process_tomography.py
index 18ae55dac8..bdc0fa3c0a 100644
--- a/test/library/tomography/test_process_tomography.py
+++ b/test/library/tomography/test_process_tomography.py
@@ -22,11 +22,11 @@
 import qiskit.quantum_info as qi
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import XGate, CXGate
-from qiskit.result import LocalReadoutMitigator
 
 from qiskit_aer import AerSimulator
 from qiskit_aer.noise import NoiseModel
 
+from qiskit_experiments.data_processing import LocalReadoutMitigator
 from qiskit_experiments.database_service import ExperimentEntryNotFound
 from qiskit_experiments.library import ProcessTomography, MitigatedProcessTomography
 from qiskit_experiments.library.tomography import ProcessTomographyAnalysis, basis
diff --git a/test/library/tomography/test_state_tomography.py b/test/library/tomography/test_state_tomography.py
index b0ec4b52c8..d4c3cb2ec4 100644
--- a/test/library/tomography/test_state_tomography.py
+++ b/test/library/tomography/test_state_tomography.py
@@ -23,11 +23,11 @@
 import qiskit.quantum_info as qi
 from qiskit import QuantumCircuit
 from qiskit.circuit.library import XGate
-from qiskit.result import LocalReadoutMitigator
 
 from qiskit_aer import AerSimulator
 from qiskit_aer.noise import NoiseModel
 
+from qiskit_experiments.data_processing import LocalReadoutMitigator
 from qiskit_experiments.database_service import ExperimentEntryNotFound
 from qiskit_experiments.library import StateTomography, MitigatedStateTomography
 from qiskit_experiments.library.tomography import StateTomographyAnalysis, basis