Fit using `BinomialFitter` failing

[AUTProgram]

I tried to use the BinomialFitter to fit SDF scans for a single ion with the MolmerSorensenFreq model. With the following code

date = "2024-01-25"
rid = 994911

data = load_result("fairy-castle", date, rid=rid)
detuning_symmetric_arr = np.array(data["datasets"][f"ndscan.rid_{rid}.points.axis_0"])
p_arr = np.array(data["datasets"][f"ndscan.rid_{rid}.points.channel_p"])
p_err_arr = np.array(data["datasets"][f"ndscan.rid_{rid}.points.channel_p_err"])

duration_estimate = 106.5e-6
n_shots = 100
model = rescale_model_x(MolmerSorensenFreq, 2 * np.pi)(num_qubits=1, start_excited=True, walsh_idx=1)
model.parameters["t_pulse"].user_estimate = duration_estimate
model.parameters["omega"].user_estimate = np.pi / duration_estimate * np.sqrt(2)
model.parameters["w_0"].fixed_to = 0.0
fit = BinomialFitter(detuning_symmetric_arr, p_arr, n_shots, model)

get this error:

C:\Users\ClemensLoschnauer\AppData\Local\pypoetry\Cache\virtualenvs\analysis-_p4OdmLl-py3.8\lib\site-packages\scipy\optimize\_numdiff.py:576: RuntimeWarning: invalid value encountered in subtract
  df = fun(x) - f0
---------------------------------------------------------------------------
RuntimeError                              Traceback (most recent call last)
Cell In[2], line 15
     13 model.parameters["omega"].user_estimate = np.pi / duration_estimate * np.sqrt(2)
     14 model.parameters["w_0"].fixed_to = 0.0
---> 15 fit = BinomialFitter(detuning_symmetric_arr, p_arr, n_shots, model)

File ~\AppData\Local\pypoetry\Cache\virtualenvs\analysis-_p4OdmLl-py3.8\lib\site-packages\ionics_fits\binomial.py:54, in BinomialFitter.__init__(self, x, y, num_trials, model, step_size)
     43 """Fits a model to a dataset and stores the results.
     44 
     45 :param x: x-axis data
   (...)
     51 :param step_size: see :class MLEFitter:
     52 """
     53 self.num_trials = num_trials
---> 54 super().__init__(x=x, y=y, model=model, step_size=step_size)

File ~\AppData\Local\pypoetry\Cache\virtualenvs\analysis-_p4OdmLl-py3.8\lib\site-packages\ionics_fits\MLE.py:60, in MLEFitter.__init__(self, x, y, model, step_size)
     57 if np.any(y < 0) or np.any(y > 1):
     58     raise RuntimeError("y values must lie between 0 and 1")
---> 60 super().__init__(x=x, y=y, model=model)

File ~\AppData\Local\pypoetry\Cache\virtualenvs\analysis-_p4OdmLl-py3.8\lib\site-packages\ionics_fits\common.py:650, in Fitter.__init__(self, x, y, model)
    647     y = scaled_model.func(x, params)
    648     return y
--> 650 fitted_params, uncertainties = self._fit(
    651     x_scaled, y_scaled, scaled_model.parameters, free_func
    652 )
    653 fitted_params.update(
    654     {param: value for param, value in self.fixed_parameters.items()}
    655 )
    656 uncertainties.update({param: 0 for param in self.fixed_parameters.keys()})

File ~\AppData\Local\pypoetry\Cache\virtualenvs\analysis-_p4OdmLl-py3.8\lib\site-packages\ionics_fits\MLE.py:130, in MLEFitter._fit(self, x, y, parameters, free_func)
    121 res = optimize.minimize(
    122     fun=self.log_liklihood,
    123     args=(x, y, free_func),
   (...)
    126     options={"maxls": 100},
    127 )
    129 if not res.success:
--> 130     raise RuntimeError(f"{self.TYPE} fit failed: {res.message}")
    132 # Compute parameter covariance matrix
    133 #
    134 # The covariance matrix is given by the inverse of the Hessian at the optimum.
   (...)
    138 # we perform our own calculation based on finite differences using the
    139 # specified step size, rescaled by the parameter bounds where possible.
    140 lower = [bound if bound is not None else -np.inf for bound in lower]

RuntimeError: Binomial fit failed: ABNORMAL_TERMINATION_IN_LNSRCH

When I use the NormalFitter with the same parameters, fit works without problems, see below.

[hartytp]

I can reproduce this with

 name         : ionics-fits
 version      : 1.3.9
 description  : Lightweight Python data fitting library with an emphasis on AMO/Quantum Information

dependencies
 - numpy >=1.21.2,<2.0.0
 - scipy >=1.7.1,<2.0.0
 - statsmodels >=0.13.2,<0.14.0

[hartytp]

After some debugging, the issue here is this: if the fit model ever tells us that the probability of a point should be 0 but the dataset has a non-zero value for that point, the log-likelihood function blows up. I'm not sure if there is a "standard" way of dealing with this, but in #177 I resolved it by replacing all points where the model function is exactly 0 or exactly 1 with a small epsilon:

eps = 1e-12
p[p == 0] = eps
p[p == 1] = 1 - eps

cc @r-srinivas