Update computation of probability

[nkanazawa1989]

Summary

As suggested by @chriseclectic and @oliverdial, new data processor node to compute probability from count dict by assuming the Dirichlet distribution is added. Conventional Probability node will be continuously supported but default node is replaced with DirichletProbability.

Details and comments

This change broke some unittests since their initial guess were really sensitive to the probability and sigma. The guess functions of DragCalAnalysis and ErrorAmplificationAnalysis have been upgraded to be more robust, and the error tolerance of the CR Hamiltonian tomography unittest has been relaxed from 1 kHz to 20 kHz (CR Hamiltonian tomography is especially sensitive to probability).

[nkanazawa1989]

It seems like the model of error amplification analysis has a strong local minima. The bug fix in 1e8750f yields the test failure (the code to take sum of alpha_i had an error), and I found that

amp_cal = FineAmplitude(qubit)
amp_cal.set_schedule(
    schedule=xp_sched, angle_per_gate=target_angle, add_xp_circuit=True, add_sx=True
)

error = -np.pi * r  # modify this line
backend = MockFineAmp(error, np.pi, "xp")

exp_data = amp_cal.run(backend)
exp_data.block_for_results()

r=0.0495 (success), 0.05 (fail), 0.0508 (success). We need almost exact d_theta estimate to avoid the minima, but this is not realistic. At least we need data driven test for fine amplitude module. Probably we need to use different fit model, or switch to global optimizer such as differential_evolution.

[oliverdial]

I'm trying to sort through the code to see what you're actually fitting. I'm not sure wrapping cos really improves readability.
https://github.com/Qiskit/qiskit-experiments/blob/1e8750f55f6645bca9733d8e4b95c34440cfe73a/qiskit_experiments/curve_analysis/fit_function.py#L21

Anyway, it looks like you're fitting both the angular error and the amplitude, but I don't see the 1/0 calibrations included in the fit. That doesn't work -- for small errors amp and dtheta_error do the same thing, so the fit function is singular. You need to either drop amp, or even better fit the 1/0 cals.

(to elucidate the later a bit; if you taylor expand for small angles, the fit function is amp * dtheta_error*(-1)^n and you can't fit both amp and detheta_error. If you have a 1/0 cal you trust (all cals except the pi/2 pulse cal), you can drop amp as a parameter, but for the pi/2 cal if the pulse is off, your calibration of what 1 and 0 are will also be off, so you're better off explicitly fitting those as well... )

[nkanazawa1989]

Thanks @oliverdial for perfectly reasonable explanation. It sounds like the problem is outside the scope of this PR, and we need to implement the error amplification family as a BatchExperiment consisting of 0/1 cal and existing error amplification experiments, and fix amp of the fit model by the value estimated by the 0/1 cal experiment. Probably a temporary solution for this PR is to set amp = 1.0 (assuming negligible assignment error).

The reason for wrapping cos function is reusability rather than readability (actually it sometimes hurts readability). If there is a set of common fit functions and parameter guess functions, we can easily propagate bug fix/upgrade through the stack.

[eggerdj]

One of the problems that we have is that the error amplifying sequences are technically not a curve analysis (i.e. not a single parameter x). Instead, it is my understanding that these error amplifying analysis classes should allow us to have a multi-dimensional x variable such as (angle_per_gate, repetition, phase_offset) that can take on different values depending on the circuit.

[eggerdj]

Would it make sense to make ErrorAmplificationAnalysis inherit from BaseAnalysis instead of CurveAnalysis? Or, would it be possible to have a curve analysis with multi-dimensional x?

[nkanazawa1989]

This means we can sweep angle_per_gate or phase_offset in a certain experiment? For me these parameters can change from experiment to experiment but what we can scan in a single experiment is only repetition. I often use fit function lambda x, b, err: ((-1) ** x * np.cos(np.pi/2 + x * err)) + b with fixed amp = 1.0 and seems like it works well. I don't think multi dimensional fit can solve the issue of simultaneous optimization for amp and d_theta (perhaps I miss something).

[oliverdial]

One approach to this is to assign arbitrary x coordinates just for the plotting part of the fitter, and then use other means to get the (richer multi-dimensional) coordinates to the fitter; I've done this a lot in the past with patterns like fit_function = lambda x, rich_x=data: actual_fit_function(rich_x)

[nkanazawa1989]

In this PR amp is fixed at 1.0 (done in d05a3b0) because this PR is originally for the introduction of new probability node. Upgrade of the fit model (to multi-dim) of error amplification will be done in a follow-up PR.

(EDIT)
This PR is ready for the final review.

[eggerdj]

Looks good. Some things could do with a bit more clarification.

[eggerdj]

I'm a bit confused by the variable name. I though mode is the value that appears most often in the data. Since we have counts I would naively expect this to be 0 or 1. Should this be named first_moment or mean instead of mode?

[nkanazawa1989]

I just referred to the definition here. According to the definition, alpha_i > 0 and it can be a negative value when alpha_i is less than 1 it says alpha_i > 1. So probably I use wrong definition. Please correct me if I miss something @oliverdial

[eggerdj]

What prompted this change?

[nkanazawa1989]

Previously, I let the FFT module estimate a frequency regardless of the signal amplitude. In some test cases, especially when the signal is really weak, the guess module picks random frequency (because computed probability has been changed), and it hurts frequency guess of averaged oscillation frequency measured in x, y, z basis. So I decided to ignore the guess when SNR is likely low, i.e. it may pick some artifact.

[oliverdial]

I like it! I wonder if using the actual SNR (min-max)/(sqrt(variance)) might be nice in case someone tries to fit a very small/off resonant rotation with this at some point, so that they can average a lot and still do it? Might be a separate PR

[nkanazawa1989]

That is nice idea. What is the variance here? Is it the variance of y values np.std(data.y) or one computed from the sampling error data.y_err? If latter one, probably np.mean(data.y_err).

[eggerdj]

Should we use the get function of the data processor library?

[nkanazawa1989]

I checked but we cannot use the getter. This is a custom processor for RB, because outcome depends on the number of qubits, and it counts for label "0" instead of "1". Of course we can modify the getter, but this should be done in the separate PR.

[eggerdj]

Can you add a few extra tests to check that we get the right probabilities? E.g. if I give it {"1": 50, "0": 50} with a flat prior I should get p = 50%.

[nkanazawa1989]

Done in e534cea This was really nice test since it told an error in the computation of mode. The bug is also fixed.

[chriseclectic]

I think this PR is confusing some things with the estimation problem and because of that the documentation and implementation is a bit confusing.

The Probability node should be estimating a single probability value p, which defines a 2-outcome binomial distribution [p, 1-p]. This is estimated from a number of samples from this distribution (the frequency of a single outcome out of total shots). The distribution over [p, 1-p] itself is a beta distribution, and so given a prior distribution of this you can do the bayes update for beta distribution to get a posterior given the observed counts, and from that distribution get a most likely value for p and a variance.

The generalization of this would be a Probabilities node which is estimating K probabilities defining a multi-normal distribution. This estimate is done from samples from this distribution (ie counts). The Dirichlet distribution is then the actual distribution over this probability vector values (generalization of beta-distribution), where you can do the Bayes update from counts the way that is coded in here and then estimate.

In this sense I think we only need two probability data processing nodes, one for a single outcome (Probability using beta distribution) and one for 2^n outcomes (Probabilities using dirichlet distribution). For the probabilities case you could return a dictionary of most likely and a separate dictionary of variances/stddev (or a single dict with pairs for values). To keep this sparse for non-observed outcomes, which for a prior of alpha=1 have MLE estimate of p_i=0 but a non-zero variance (and a non-zero est for p_i for others priors) you could have some sort of default value dict (i think there is something like this in python collections?).

Now if the scope of this PR is just to be updating the Probability node to compute standard deviation correctly then you can just update that node to have corrected API docs referring to beta distribution and update the calculation accordingly. The full Probabilities node could be added later if needed (Tomography in principle would use this data processor but at the moment has its own version without using the data processor classes)

[chriseclectic]

I think these docs are getting a little too detailed with examples for API docs and getting more into user guide territory, which we should eventually be separating into its own section.

[nkanazawa1989]

I think these description is necessary to explain what is alpha_i. If you think saying alpha_i is a parameter of Dirichlet distribution is enough, probably these details can be removed (but anyways I don't think there is a place to write userguide of dirichlet distribution since this is not an experiment).

[chriseclectic]

After getting to here I think the whole distribution stuff in the into is a bit confusion now. The Dirichlet distribution if for computing a multi-variate probability distribution, so for example converting observed counts into an estimate of a probability distribution over all possible outcomes (2^n). If it is for a binary estimate this should be a beta distribution (which is the special case of 2-outcome dirichlet distribution).

The way the description of this class is worded it should be something that converts a full counts dict into a full probabilities dict + std error dict since the bayesian update is for a full distribution on all outcomes, not a distribution for a single outcome. The difference between estimating individual probabilities as 2-outcome distributions vs all probabilities as 2^n outcome distribution will result in slight differences in variances and what they mean.

[nkanazawa1989]

I agree the documentation for the option is confusing, i.e. it uses multivariable distribution but distribution defaults to 2-outcome. I updated the documentation and removed the default value in 194568b

Addressed