def getPrimes(n):
    primes = list(range(2, n + 1))
    for i in primes:
        if i:
            for j in range(2 * i - 2, n - 1, i):
                primes[j] = 0

    return [i for i in primes if i]

def probability(n, p, k):
    return sps.binom(n, k) * (p ** k) * ((1 - p) ** (n - k))

def f(n, p):
    return sum([prob(n, p, k) for k in getPrimes(n)])


getPrimes считает все простые от 2 до n включительно
probability считает вероятность P(X = k) X~Binomial(n, p)
f считает сумму по всем простым этих вероятностей

for n in range(10, 101, 10):
    prs = [round(primes_p(n, p), 3) for p in np.linspace(0.1, 0.9, 9)]
    print(prs)
этим циклом можно посчитать и вывести вероятности из f так, что по вдоль строки p меняется от 0.1 до 0.9 а вдоль столбцов n меняется от 10 до 100 с шагом в 10

fig = plt.figure(1, figsize=(10, 5))
axes = fig.add_subplot(1, 1, 1)

[plt.plot(
    range(10, 101, 10),
    [round(primes_p(n, p), 3) for n in range(10, 101, 10)],
    label = 'p = ' + str(round(p, 1))
) for p in np.linspace(0.1, 0.9, 9)]

axes.legend()
plt.show()

а вот так можно все эти штуки нарисовать