# max_shots and max_errors to reduce uncerainty in threshold curves with sinter / stim?

thanks for all your answers so far! I'm wondering, what do max_shots and max_errors mean (as in the code below) when sampling from a stabiliser circuit to produce threshold curves using stim? Would changing their values reduce uncertainty when plotting threshold curves for quite large distances of the surface code? The code I'm using to produce the graph is from Stim's Getting Started notebook but modified to include higher distances:

import scipy.stats
import numpy as np
import stim
import pymatching
import sinter
from typing import *
import matplotlib.pyplot as plt

circuit = stim.Circuit.generated(
"surface_code:rotated_memory_z",
rounds=d * 3,
distance=d,
after_clifford_depolarization=noise,
after_reset_flip_probability=noise,
before_measure_flip_probability=noise,
before_round_data_depolarization=noise,
),
json_metadata={'d': d, 'r': d * 3, 'p': noise},
)
for d in [3, 5, 7, 9, 11, 13,15,17,19,21]
for noise in [0.008, 0.009, 0.01, 0.011, 0.012,0.013]
]

num_workers=4,
decoders=['pymatching'],
max_shots=10_000_000,
max_errors=5_000,
print_progress=False,
)

#Plot logical error rate per round of stab measurements for each distance

fig, ax = plt.subplots(1, 1)
sinter.plot_error_rate(
ax=ax,
stats=collected_surface_code_stats,
)
# ax.set_ylim(5e-3, 5e-2)
# ax.set_xlim(0.008, 0.02)
ax.loglog()
ax.set_title("Rotated Surface Code Error Rates per Round under Circuit Noise")
ax.set_xlabel("Physical Error Rate")
ax.set_ylabel("Logical Error Rate per Round")
ax.grid(which='major')
ax.grid(which='minor')
ax.legend()
fig.set_dpi(120)  # Show it bigger


The underlying issue here is that the meaning of "per-round error" changes discontinuously as the hypothetical per-shot error rate crosses 50%, and you are looking at samples where the uncertainty in the per-shot error rate crosses 50%. The argument failure_units_per_shot_func=lambda stats: stats.json_metadata['r'] is telling sinter to convert from per-shot to per-round error rates and display those. It is extremely difficult to fix this issue by taking more shots; it's sort of intrinsic to doing the conversion.

For example, suppose the experiment had 1000 rounds and you took two million shots and you saw exactly one million errors. The sample mean is 50%. The sample standard deviation is ~0.03%, so let's say the uncertainty region is 49.7% to 50.3%. An error rate above 50% sounds a bit odd, and as we will see it also behaves a bit oddly, but in principle it can actually happen. For example, the decoder might be accidentally inverting its answers.

Because there are 1000 rounds, the per-round error rate $$p_\text{round}$$ must satisfy the property that, if you sampled a Bernoulli distribution 1000 times, the probability of there being an odd number of errors is $$p_\text{shot}$$. For the lower bound $$p_\text{shot}=49.7\%$$ you can calculate the corresponding $$p_\text{round} \approx 0.25\%$$. For the upper bound $$p_\text{shot}=50.3\%$$ you can calculate the corresponding $$p_\text{round} \approx 99.75\%$$.

See the problem? Even though we took a lot of shots and saw an absolutely enormous number of errors and successes, the conversion is amplifying the small per-shot uncertainty into an enormous per-round uncertainty. A 0.06% gap became a 99.5% gap! Basically, we're unable to tell whether or not the decoder is accidentally inverting results and this completely confounds what "per round error rate" means.

Sinter could assume that error rates were promised to be below 50% in order to avoid showing this gap, or define per-round error rates to always be below 50%. I'm not sure I feel comfortable with making those kinds of assumptions automatically. You can in fact create a round that has a greater than 50% error rate, by inserting a logical bit flip 100% of the time. Showing a giant gross region clearly indicating the statistical difficulty is the current compromise.

As an aside, I'm personally of the opinion that the subtleties of the behavior of logical errors at high physical error rates don't matter. They're irrelevant, if not actively misleading, to predicting costs in the actually-working regime. For a 1000 round experiment, per-round error rates of 1%, 10%, 50%, and 90% all correspond to same thing: "the computation basically never succeeds". You can't build a computer that way. I actually kinda like what the diagram is doing. Not wanting that part of the plot to look like that is a great metaphor for not wanting the physical noise rate to be that bad.

Appendix: combining Bernoulli distributions.

To xor the outputs of two Bernoulli distributions means to sample them both and return whether they are the same (OFF) or different (ON). If the two distributions have probability-of-ON $$p_1$$ and $$p_2$$, then their xor is also a Bernoulli distribution and it has probability-of-ON $$p_{1 \oplus 2} = p_1 \overline{p_2} + p_2 \overline{p_1}$$ where $$\overline{x} = 1 - x$$. When $$p_1=p_2$$ this simplifies to $$p_{\oplus} = 2 p \overline{p} = 2p - 2p^2$$.

A particularly useful conversion here is to switch from probability-of-$$ON$$ $$p$$ to "probability of preservation" $$t$$ where $$t=\overline{2p}$$. In this view, the Bernoulli distribution is a communication channel and there is a probability $$t$$ of the message being transmitted successfully and a probability $$\overline{t}$$ of it being replaced by a random bit. The self-xoring equation is much simpler from this perspective: $$t_{\oplus} = t^2$$. In fact, this formula generalizes. If you xor a distribution with itself $$r$$ times, the preservation rate of the resulting distribution is $$t^r$$.

So, we're converting from per-shot error rates $$p_s$$ to per-round error rates $$p_r$$. The definition being used is that the per-shot error rate is what you get from xoring $$r$$ per-round Bernoulli distributions together. By converting to preservation probabilities and figuring out the relationship there then converting back, you get the equation

$$p_r = \frac{1 - (1 - 2p_s)^{1/r}}{2}$$

To extend this equation to $$p_s > 50\%$$, repeat the whole probability-of-preservation logic but with probability-of-sending-flipped-message. You then get:

$$p_r = f_r(p_s) = \begin{cases}p_s > 50\% & 1 - f_r(1-p_s)\\p_s \leq 50\% & \frac{1 - (1 - 2p_s)^{1/r}}{2}\end{cases}$$