Update computation of probability

qiskit_experiments/curve_analysis/standard_analysis/error_amplification_analysis.py

@@ -17,7 +17,6 @@
 import numpy as np
 
 import qiskit_experiments.curve_analysis as curve
-from qiskit_experiments.exceptions import CalibrationError
 
 
 class ErrorAmplificationAnalysis(curve.CurveAnalysis):
@@ -32,7 +31,29 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
 
         .. math::
             y = \frac{{\rm amp}}{2}\cos\left(x[{\rm d}\theta + {\rm apg} ] \
-            +{\rm phase\_offset}\right)+{\rm base}
+            -{\rm phase\_offset}\right)+{\rm base}
+
+        To understand how the error is measured we can transformed the function above into
+
+        .. math::
+            y = \frac{{\rm amp}}{2} \left(\
+            \cos\right({\rm d}\theta \cdot x\left)\
+            \cos\right({\rm apg} \cdot x - {\rm phase\_offset}\left) -\
+            \sin\right({\rm d}\theta \cdot x\left)\
+            \sin\right({\rm apg} \cdot x - {\rm phase\_offset}\left)
+            \right) + {\rm base}
+
+        When :math:`{\rm apg} \cdot x - {\rm phase\_offset} = (2n + 1) \pi/2` is satisfied the
+        fit model above simplifies to
+
+        .. math::
+            y = \mp \frac{{\rm amp}}{2} \sin\left({\rm d}\theta \cdot x\right) + {\rm base}
+
+        In the limit :math:`{\rm d}\theta \ll 1`, the error can be estimated from the curve data
+
+        .. math::
+            {\rm d}\theta \simeq \mp \frac{2(y - {\rm base})}{x \cdot {\rm amp}}
+
 
     # section: fit_parameters
         defpar \rm amp:
@@ -49,6 +70,8 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
             desc: The angle offset in the gate that we wish to measure.
             init_guess: Multiple initial guesses are tried ranging from -a to a
                 where a is given by :code:`max(abs(angle_per_gate), np.pi / 2)`.
+                Extra guesses are added based on curve data when either :math:`\rm amp` or
+                :math:`\rm base` is :math:`\pi/2`. See fit model for details.
             bounds: [-pi, pi].
 
     # section: note
@@ -63,7 +86,7 @@ class ErrorAmplificationAnalysis(curve.CurveAnalysis):
                 x,
                 amp=0.5 * amp,
                 freq=(d_theta + angle_per_gate) / (2 * np.pi),
-                phase=phase_offset,
+                phase=-phase_offset,
                 baseline=base,
             ),
             plot_color="blue",
@@ -86,7 +109,6 @@ def _default_options(cls):
                 :math:`\pi/2` if the square-root of X gate is added before the repeated gates.
                 This is decided for the user in :meth:`set_schedule` depending on whether the
                 sx gate is included in the experiment.
-            number_of_guesses (int): The number of initial guesses to try.
             max_good_angle_error (float): The maximum angle error for which the fit is
                 considered as good. Defaults to :math:`\pi/2`.
         """
@@ -96,7 +118,6 @@ def _default_options(cls):
         default_options.ylabel = "Population"
         default_options.angle_per_gate = None
         default_options.phase_offset = 0.0
-        default_options.number_guesses = 21
         default_options.max_good_angle_error = np.pi / 2
 
         return default_options
@@ -115,8 +136,6 @@ def _generate_fit_guesses(
         Raises:
             CalibrationError: When ``angle_per_gate`` is missing.
         """
-        n_guesses = self._get_option("number_guesses")
-
         curve_data = self._data()
         max_abs_y, _ = curve.guess.max_height(curve_data.y, absolute=True)
         max_y, min_y = np.max(curve_data.y), np.min(curve_data.y)
@@ -128,15 +147,35 @@ def _generate_fit_guesses(
             user_opt.p0.set_if_empty(amp=max_y - min_y)
             user_opt.bounds.set_if_empty(amp=(-2 * max_abs_y, 2 * max_abs_y))
 
-        # Base the initial guess on the intended angle_per_gate.
-        angle_per_gate = self._get_option("angle_per_gate")
-
-        if angle_per_gate is None:
-            raise CalibrationError("The angle_per_gate was not specified in the analysis options.")
+        # Base the initial guess on the intended angle_per_gate and phase offset.
+        apg = self._get_option("angle_per_gate")
+        phi = self._get_option("phase_offset")
+
+        # Prepare logical guess for specific condition (often satisfied)
+        d_theta_guesses = []
+
+        offsets = apg * curve_data.x + phi
+        amp = user_opt.p0.get("amp", self._get_option("amp"))
+        for i in range(curve_data.x.size):
+            xi = curve_data.x[i]
+            yi = curve_data.y[i]
+            if np.isclose(offsets[i] % np.pi, np.pi / 2) and xi > 0:
+                # Condition satisfied: i.e. cos(apg x - phi) = 0
+                err = -np.sign(np.sin(offsets[i])) * (yi - user_opt.p0["base"]) / (0.5 * amp)
+                # Validate estimate. This is just the first order term of Maclaurin expansion.
+                if np.abs(err) < 0.5:
+                    d_theta_guesses.append(err / xi)
+                else:
+                    # Terminate guess generation because larger d_theta x will start to
+                    # reduce net y value and underestimate the rotation.
+                    break
+
+        # Add naive guess for more coverage
+        guess_range = max(abs(apg), np.pi / 2)
+        d_theta_guesses.extend(np.linspace(-guess_range, guess_range, 11))
 
-        guess_range = max(abs(angle_per_gate), np.pi / 2)
         options = []
-        for d_theta_guess in np.linspace(-guess_range, guess_range, n_guesses):
+        for d_theta_guess in d_theta_guesses:
             new_opt = user_opt.copy()
             new_opt.p0.set_if_empty(d_theta=d_theta_guess)
             options.append(new_opt)

qiskit_experiments/data_processing/nodes.py

@@ -13,7 +13,8 @@
 """Different data analysis steps."""
 
 from abc import abstractmethod
-from typing import Any, Dict, List, Optional, Tuple, Union
+from numbers import Number
+from typing import Any, Dict, List, Optional, Tuple, Union, Sequence
 import numpy as np
 
 from qiskit_experiments.data_processing.data_action import DataAction, TrainableDataAction
@@ -391,17 +392,72 @@ def _process(self, datum: np.array, error: Optional[np.array] = None) -> np.arra
 
 
 class Probability(DataAction):
-    """Count data post processing. This returns the probabilities of the outcome string
-    used to initialize an instance of Probability."""
+    r"""Compute the mean probability of a single measurement outcome from counts.
 
-    def __init__(self, outcome: str = "1", validate: bool = True):
+    This node returns the mean and standard deviation of a single measurement
+    outcome probability $p$ estimated from the observed counts. The mean and
+    variance are computed from the posterior Beta distribution
+    $B(\alpha_0^\prime,\alpha_1^\prime)$ estimated from a Bayesian update
+    of a prior Beta distribution $B(\alpha_0, \alpha_1)$ given the observed
+    counts.
+
+    The mean and variance of the Beta distribution $B(\alpha_0, \alpha_1)$ are:
+
+    .. math::
+
+        \text{E}[p] = \frac{\alpha_0}{\alpha_0 + \alpha_1}, \quad
+        \text{Var}[p] = \frac{\text{E}[p] (1 - \text{E}[p])}{\alpha_0 + \alpha_1 + 1}
+
+    Given a prior Beta distribution $B(\alpha_0, \alpha_1)$, the posterior
+    distribution for the observation of $F$ counts of a given
+    outcome out of $N$ total shots is a $B(\alpha_0^\prime,\alpha_1^\prime)$ with
+
+    .. math::
+        \alpha_0^\prime = \alpha_0 + F, \quad
+        \alpha_1^\prime = \alpha_1 + N - F.
+
+    .. note::
+
+        The default value for the prior distribution is *Jeffery's Prior*
+        $\alpha_0 = \alpha_1 = 0.5$ which represents ignorance about the true
+        probability value. Note that for this prior the mean probability estimate
+        from a finite number of counts can never be exactly 0 or 1. The estimated
+        mean and variance are given by
+
+        .. math::
+
+            \text{E}[p] = \frac{F + 0.5}{N + 1}, \quad
+            \text{Var}[p] = \frac{\text{E}[p] (1 - \text{E}[p])}{N + 2}
+    """
+
+    def __init__(
+        self,
+        outcome: str,
+        alpha_prior: Union[float, Sequence[float]] = 0.5,
+        validate: bool = True,
+    ):
         """Initialize a counts to probability data conversion.
 
         Args:
-            outcome: The bitstring for which to compute the probability which defaults to "1".
+            outcome: The bitstring for which to return the probability and variance.
+            alpha_prior: A prior Beta distribution parameter ``[`alpha0, alpha1]``.
+                         If specified as float this will use the same value for
+                         ``alpha0`` and``alpha1`` (Default: 0.5).
             validate: If set to False the DataAction will not validate its input.
+
+        Raises:
+            DataProcessorError: When the dimension of the prior and expected parameter vector
+                do not match.
         """
         self._outcome = outcome
+        if isinstance(alpha_prior, Number):
+            self._alpha_prior = [alpha_prior, alpha_prior]
+        else:
+            if validate and len(alpha_prior) != 2:
+                raise DataProcessorError(
+                    "Prior for probability node must be a float or pair of floats."
+                )
+            self._alpha_prior = list(alpha_prior)
         super().__init__(validate)
 
     def _format_data(self, datum: dict, error: Optional[Any] = None) -> Tuple[dict, Any]:
@@ -472,12 +528,14 @@ def _process(
 
             return np.array(populations), np.array(errors)
 
-    def _population_error(self, counts_dict) -> Tuple[float, float]:
+    def _population_error(self, counts_dict: Dict[str, int]) -> Tuple[float, float]:
         """Helper method"""
         shots = sum(counts_dict.values())
-        p_mean = counts_dict.get(self._outcome, 0.0) / shots
-        p_var = p_mean * (1 - p_mean) / shots
-
+        freq = counts_dict.get(self._outcome, 0)
+        alpha_posterior = [freq + self._alpha_prior[0], shots - freq + self._alpha_prior[1]]
+        alpha_sum = sum(alpha_posterior)
+        p_mean = alpha_posterior[0] / alpha_sum
+        p_var = p_mean * (1 - p_mean) / (alpha_sum + 1)
         return p_mean, np.sqrt(p_var)
 
 

qiskit_experiments/data_processing/processor_library.py

@@ -16,7 +16,7 @@
 
 from qiskit_experiments.data_processing.exceptions import DataProcessorError
 from qiskit_experiments.data_processing.data_processor import DataProcessor
-from qiskit_experiments.data_processing.nodes import AverageData, Probability, SVD, MinMaxNormalize
+from qiskit_experiments.data_processing import nodes
 
 
 def get_processor(
@@ -38,16 +38,16 @@ def get_processor(
         DataProcessorError: if the measurement level is not supported.
     """
     if meas_level == MeasLevel.CLASSIFIED:
-        return DataProcessor("counts", [Probability("1")])
+        return DataProcessor("counts", [nodes.Probability("1")])
 
     if meas_level == MeasLevel.KERNELED:
         if meas_return == "single":
-            processor = DataProcessor("memory", [AverageData(axis=1), SVD()])
+            processor = DataProcessor("memory", [nodes.AverageData(axis=1), nodes.SVD()])
         else:
-            processor = DataProcessor("memory", [SVD()])
+            processor = DataProcessor("memory", [nodes.SVD()])
 
         if normalize:
-            processor.append(MinMaxNormalize())
+            processor.append(nodes.MinMaxNormalize())
 
         return processor
 

qiskit_experiments/library/calibration/analysis/drag_analysis.py

@@ -156,16 +156,19 @@ def _generate_fit_guesses(
             beta=(-freq_bound, freq_bound),
             base=(-max_abs_y, max_abs_y),
         )
-        user_opt.p0.set_if_empty(amp=-0.5, base=0.5)
+        user_opt.p0.set_if_empty(base=0.5)
 
         # Drag curves can sometimes be very flat, i.e. averages of y-data
         # and min-max do not always make good initial guesses. We therefore add
-        # 0.5 to the initial guesses.
+        # 0.5 to the initial guesses. Note that we also set amp=-0.5 because the cosine function
+        # becomes +1 at zero phase, i.e. optimal beta, in which y data should become zero
+        # in discriminated measurement level.
         options = []
-        for beta_guess in np.linspace(min_beta, max_beta, 20):
-            new_opt = user_opt.copy()
-            new_opt.p0.set_if_empty(beta=beta_guess)
-            options.append(new_opt)
+        for amp_guess in (0.5, -0.5):
+            for beta_guess in np.linspace(min_beta, max_beta, 20):
+                new_opt = user_opt.copy()
+                new_opt.p0.set_if_empty(amp=amp_guess, beta=beta_guess)
+                options.append(new_opt)
 
         return options
 

qiskit_experiments/library/calibration/analysis/fine_amplitude_analysis.py

@@ -28,4 +28,6 @@ class FineAmplitudeAnalysis(ErrorAmplificationAnalysis):
     """
 
     # The intended angle per gat of the gate being calibrated, e.g. pi for a pi-pulse.
-    __fixed_parameters__ = ["angle_per_gate", "phase_offset"]
+
+    # TODO remove amp from fixed parameter.
+    __fixed_parameters__ = ["angle_per_gate", "phase_offset", "amp"]

qiskit_experiments/library/calibration/fine_drag.py

@@ -160,7 +160,7 @@ def _default_analysis_options(cls) -> Options:
         options = super()._default_analysis_options()
         options.normalization = True
         options.angle_per_gate = 0.0
-        options.phase_offset = -np.pi / 2
+        options.phase_offset = np.pi / 2
         options.amp = 1.0
 
         return options

qiskit_experiments/library/characterization/cr_hamiltonian_analysis.py

@@ -199,7 +199,8 @@ def _default_options(cls):
         """Return the default analysis options."""
         default_options = super()._default_options()
         default_options.data_processor = dp.DataProcessor(
-            input_key="counts", data_actions=[dp.Probability("1"), dp.BasisExpectationValue()]
+            input_key="counts",
+            data_actions=[dp.Probability("1"), dp.BasisExpectationValue()],
         )
         default_options.curve_plotter = "mpl_multiv_canvas"
         default_options.xlabel = "Flat top width"
@@ -266,21 +267,23 @@ def _generate_fit_guesses(
 
         guesses = defaultdict(list)
         for control in (0, 1):
-            # start from Z oscillation
             x_data = self._data(series_name=f"x|c={control}")
             y_data = self._data(series_name=f"y|c={control}")
             z_data = self._data(series_name=f"z|c={control}")
 
             omega_xyz = []
             for data in (x_data, y_data, z_data):
+                ymin, ymax = np.percentile(data.y, [10, 90])
+                if ymax - ymin < 0.2:
+                    # oscillation amplitude might be almost zero,
+                    # then exclude from average because of lower SNR
+                    continue
                 fft_freq = curve.guess.frequency(data.x, data.y)
-                # oscillation amplitude might be almost zero, then exclude from average
-                if fft_freq > 0:
-                    omega_xyz.append(fft_freq)
+                omega_xyz.append(fft_freq)
             if omega_xyz:
                 omega = 2 * np.pi * np.average(omega_xyz)
             else:
-                omega = 0.0
+                omega = 1e-3
 
             zmin, zmax = np.percentile(z_data.y, [10, 90])
             theta = np.arccos(np.sqrt((zmax - zmin) / 2))

qiskit_experiments/library/characterization/fine_amplitude.py

@@ -118,6 +118,14 @@ def _default_experiment_options(cls) -> Options:
 
         return options
 
+    @classmethod
+    def _default_analysis_options(cls) -> Options:
+        """Default analysis options."""
+        options = super()._default_analysis_options()
+        options.amp = 1.0
+
+        return options
+
     def __init__(self, qubit: int, gate: Gate):
         """Setup a fine amplitude experiment on the given qubit.
 

qiskit_experiments/library/randomized_benchmarking/rb_experiment.py

@@ -22,8 +22,8 @@
 from qiskit.quantum_info import Clifford
 from qiskit.circuit import Gate
 
+import qiskit_experiments.data_processing as dp
 from qiskit_experiments.framework import BaseExperiment, ParallelExperiment, Options
-from qiskit_experiments.curve_analysis.data_processing import probability
 from .rb_analysis import RBAnalysis
 from .clifford_utils import CliffordUtils
 from .rb_utils import RBUtils
@@ -88,7 +88,12 @@ def __init__(
 
         # Set configurable options
         self.set_experiment_options(lengths=list(lengths), num_samples=num_samples)
-        self.set_analysis_options(data_processor=probability(outcome="0" * self.num_qubits))
+        self.set_analysis_options(
+            data_processor=dp.DataProcessor(
+                input_key="counts",
+                data_actions=[dp.Probability(outcome="0" * self.num_qubits)],
+            )
+        )
 
         # Set fixed options
         self._full_sampling = full_sampling

releasenotes/notes/upgrade-probability-computation-6c7539a52e3d4484.yaml

@@ -0,0 +1,14 @@
+---
+fixes:
+  - |
+    Update the :class:`~qiskit_experiments.data_processing.Probability`
+    data processing node to compute the estimated mean and standard deviation
+    of a measurement outcome probability using a Bayesian update of a
+    a Beta distribution prior given the observed measurement outcomes.
+
+    The default uninformative prior assumes ignorance about the probability
+    to be estimated and will prevent the estimated mean from being exactly
+    0 or 1, and prevent the estimated standard deviation from being 0, which
+    could cause issues with computing weights during curve fitting. A custom
+    prior can also be provided if prior information about the probability is
+    know.