The reason why Agresti and Coull chose these "add two successes and two failures" lies in rounding of a Wilson 95% CI. It was not determined by simulation studies.
I had to write a small paper on the topic once, here were my findings (condensed) and an answer to your question.
Pro-tip, R has a package propCIs which has different CI's build in. (Agresti-Coull is the add4ci method)
Background
When $X_i \stackrel{d}{=} \text{Ber}(p)$ and the sample size is $n$ then one would intuitively try the following $1-\alpha$ CI (so called Wald CI):
$$\hat p \pm z_{\alpha/2}\sqrt{\frac{\hat p(1-\hat p)}{n}}$$
But as you know this behaves badly.
One of the reason the default CI behaves badly, is the usage of $\hat p$ to find the width of the CI. This results in a very small width when $\hat p$ is close to 1 or 0. The Wilson CI tries to mitigate this issue.
This Wilson CI uses the width under $H_0$ and seeks all $p$ which solve
$$\left|\dfrac{\hat p - p}{\sqrt{\frac{p(1-p)}{n}}}\right| < z_{\alpha/2}$$
Working this out results in a quadratic equation wich results in the following crazy formula for the CI:
$$\hat p \left(\frac{n}{n+z_{\alpha/2}^2}\right) + \frac{1}{2}\left( \frac{z_{\alpha/2}^2}{n+z_{\alpha/2}^2}\right) \pm z_{\alpha/2} \sqrt{\dfrac{1}{n+z_{\alpha/2}^2}\left[ \hat p (1-\hat p) \left( \dfrac{n}{n+z_{\alpha/2}^2}\right)+\dfrac{1}{2}\left(1-\dfrac{1}{2}\right) \left(\dfrac{z_{\alpha/2}^2}{n+z_{\alpha/2}^2}\right)\right]}.$$
This CI behaves pretty good. See the figure below.
Agresti-Coull
Agresti and Coull looked at the center of the Wilson CI and noticed a simplification if one calculates a 95% CI. $z_{0.025} = 1.96\approx 2$.
Now notice how the center of the Wilson CI was given by:
$$\hat p \left( \dfrac{n}{n+z^2_{\alpha/2}}\right) + \dfrac{1}{2}\left( \dfrac{z^2_{\alpha/2}}{n+z^2_{\alpha/2}} \right)$$
When you apply the simplification suggested above you find:
$$\tilde p = \hat p \left( \frac{n}{n+4}\right) + \dfrac{1}{2}\left( \dfrac{4}{n+4}\right) = \hat p \left( \frac{n}{n+4}\right) + \dfrac{2}{n+4} = \dfrac{X+2}{n+4}$$
Which explaines the "adding two failures, two success - method".
The Agresti-Coull CI is then defined as:
$$\tilde p \pm z_{\alpha/2}\sqrt{\dfrac{\tilde p(1-\tilde p)}{\tilde n}}$$
Comparison - Coverage probability
The following picture show the three (as well as Clopper-Pearson) and the coverage probability for simultation of 5000 times.
Why not use 3 or even more successes/failures ratio's?
First of all, the derivation of the Agresti-Coull interval makes sense.
I've looked at your graph of the performance of the different methods and you claim the "adding three successes" is better, but I'm not convinced. I would say it is worse since the coverage probability is systematically to large. Meaning that the CI are to large, which makes them to conservative. Using this interval it would be harder to detect a significant result.
The paper L. D. Brown et al. (2001). Interval Estimation for a Binomial Proportion. Statistical Science. contains a very good overview of all the CI's for a binomial distribution.
