import numpy as np
import scipy.stats as st
from scipy.optimize import minimize
import matplotlib.pyplot as plt

N = 101
K = 10
p = 0.75
B = st.binom(K, p)
np.random.seed(291)
x = B.rvs(size = N)

def ll_binomial(theta, data):
    sx = np.sum(data["x"])
    K = data["K"]
    N = data["N"]
    l1p = np.log1p(-theta) # np.log(1 + -p)
    out = np.log(theta) * sx + l1p * (N * K - sx)
    return -out

rng = np.random.default_rng()
init = rng.random()
data = {"x": x, "K": K, "N": np.size(x)}
minimize(ll_binomial, init, args = (data,), bounds = [(1e-5, 1 - 1e-5)])

  message: CONVERGENCE: NORM OF PROJECTED GRADIENT <= PGTOL
  success: True
   status: 0
      fun: 576.050536971854
        x: [ 7.426e-01]
      nit: 7
      jac: [ 0.000e+00]
     nfev: 18
     njev: 9
 hess_inv: <1x1 LbfgsInvHessProduct with dtype=float64>