confidence interval for biased coin's success rate depending on number of trials I have a biased coin where success is defined as head and the (un)known success rate is 0.75. Some clever mathematician told me about a formula:

which I replicated in the R code below (it should produce the CI based on alpha 0.05 I think). I also include the graph. The question I have, can this be correct as the CI goes above 1 and below 0 for let us say a sample size of 1? Maybe my code is wrong and/or the formula? Of course I hard coded the success rate and one cannot achieve a p hat of 0.75 for 1 or 2 trials. In the real world this would fluctuate. Thanks!

library(dplyr)
library(ggplot2)

df <- data.frame(
    number_of_samples = seq(1, 100, by=1)
    , success_rate = 0.75
) %>%
mutate(
    upper_ci = success_rate + (1.96 * sqrt((success_rate * (1 - success_rate))/number_of_samples))
    , lower_ci = success_rate - (1.96 * sqrt((success_rate * (1 - success_rate))/number_of_samples))
)

font_size = 12
options(repr.plot.width = 7, repr.plot.height = 5)
ggplot() + 
    ggtitle("") +
    xlab("") +
    ylab("") +
    geom_line(data=df, aes(x=number_of_samples, y=success_rate), size=1, color = "red") +
    geom_line(data=df, aes(x=number_of_samples, y=upper_ci), size=1, color = "blue") +
    geom_line(data=df, aes(x=number_of_samples, y=lower_ci), size=1, color = "blue") +
    theme(
        axis.title.y=element_blank(),
        #axis.text.y = element_text(color="black", size = 10, angle = 0, hjust = .5, vjust = .5, face = "bold"),
        axis.text.x=element_text(angle = 0, color="black", size = font_size, face = "bold", hjust = .5, vjust = .5),
        axis.text.y=element_text(color="black", size = font_size, face = "bold", hjust = .5, vjust = .5),
        #aspect.ratio=1,
        legend.title=element_blank(),
        text = element_text(size=font_size, face="bold"),
        plot.title = element_text(hjust = 0.5)
    ) +
    coord_cartesian(ylim = c(0, 2)) +
    scale_y_continuous(breaks = round(seq(0, 2, by = 0.1), 1))

 A: As the OP correctly observed, the CI for $\hat{p}$ should not go beyond zero or one, as it is an estimate for the probability of success.
The main reason that the CI in the question beaches these two thresholds is because it uses a formula that assumes $\hat{p}$ is (approximately) normally distributed. While the approximation is sound for large $n$ (and/or $p$ close to 0.5), it does not hold for small $n$.
To understand the behaviour of $\hat{p}$ and its associated CI, we have to look at the actual, underlying generating process. If we flip a coin with success rate $\hat{p}$ for $n$ times, then the number of successes follows a binomial distribution $Bin(n, \hat{p})$. We can then obtain the CI of $\hat{p}$ as the $\alpha/2$ and $1-\alpha/2$ quantile of the said binomial distribution.
Python code:
import numpy as np
from scipy.stats import binom
from matplotlib import pyplot as plt

p = 0.75
alpha = 0.05

n = np.arange(1, 101)
CI_high = binom.ppf(1 - alpha/2, n, p)
CI_low = binom.ppf(alpha/2, n, p)

plt.plot(n, CI_high, c='blue')
plt.plot(n, CI_low, c='blue')
plt.plot(n, n*p, c='red')
plt.xlabel("Trials (n)")
plt.ylabel("Number of successes")

plt.show()

R code:
library(ggplot2)

p <- 0.75
alpha <- 0.05

n <- seq(1, 100, by=1)

CI_high <- qbinom(p=1-alpha/2, size=n, prob=p)
CI_low <- qbinom(p=alpha/2, size=n, prob=p)

df_trials <- data.frame(n, CI_high, CI_low)

ggplot() +
    xlab("Trials (n)") + ylab("Number of successes") +
    geom_line(data=df_trials, aes(x=n, y=CI_high), color='blue') +
    geom_line(data=df_trials, aes(x=n, y=CI_low), color='blue') +
    geom_line(data=df_trials, aes(x=n, y=n*p), color='red')

R plot (Blue lines are the CI, Red line is expected number of successes, i.e. $n \hat{p}$):

All three lines are increasing (at slightly different rate), which is expected. Of course, we are interested in the success rate, which plot can be easily obtained by dividing the lines by $n$:
Python code (continuation of Python code above):
plt.plot(n, CI_high/n, c='blue')
plt.plot(n, CI_low/n, c='blue')
plt.plot(n, [p] * len(n), c='red')
plt.xlabel("Trials (n)")
plt.ylabel("Success rate (p)")

plt.show()

R code (Continuation of R code above):
df_success_rate <- data.frame(n, CI_high = CI_high / n, CI_low = CI_low / n)

ggplot() +
    xlab("Trials (n)") + ylab("Success rate (p)") +
    geom_line(data=df_success_rate, aes(x=n, y=CI_high), color='blue') +
    geom_line(data=df_success_rate, aes(x=n, y=CI_low), color='blue') +
    geom_line(data=df_success_rate, aes(x=n, y=p), color='red')

...which gives the following R plot:

One can observe the shape of the CI curves are quite similar to that in the original question, except it is no longer going beyond zero or one. Note as binomials are a discrete r.v., to guarantee a $1-\alpha$ CI coverage you sometimes need the same number of successful trials for certain $n$ as compared to that for $n-1$, this leads to the rugged shape of the CI curves, as the denominator is different.
