In this answer there is a simple method given to compute the "worst case pseudo-threshold" of the Shor code. I want to apply this to the $[[5,1,3]]$ perfect code. Namely, given a depolarizing channel (where the probability that $X$, $Y$ ,or $Z$ happens is $p/3$ and the probability that $I$ happens is $1-p$) we can compute the probability of a logical error as $$ p_L^\text{worst-case} = 1-( (1-p)^5 + 5 p (1-p)^4). $$ Basically the reasoning is that because the code has distance $d=3$ we can correct weight-0 errors (which happen with probability $(1-p)^5$) and we can correct weight-1 errors (which happen with probability $p(1-p)^4$ and there are 5 of them). Then we take the complement to find the probability of a logical error.
However, this is a "worst-case" calculation because it assumes anytime we get an error above weight-1 the code fails, which isn't necessarily true (as mentioned in this answer to my previous question).
On the other hand (and mentioned in the previous question) we can try to find $p_L$ by simulating the decoding of the $[[5,1,3]]$ code in Mathematica (I used 30000 trials because it seemed smooth enough). I have plotted the probability of a logical error (y-axis) vs probability of a physical error $p$ (x-axis) for both $p_L^\text{worst-case}$ (red) and $p_L^\text{simulation}$ (blue) below. We clearly see that $p_L^\text{worst-case}$ approximates $p_L^\text{simulation}$ pretty well for small values of $p$ but deviates significantly for higher values of $p$.

I am interested in analytically deriving $p_L$ (the blue line) from first principles.
