I am estimating an Interval Regression model. I wrote my likelihood function (see intreg, STATA) as stated bellow, and compared my results with STATA. I am finding almost exactly (within 0.001) same parameters as STATA but the BFGS algorithm seems still not to converge.

I still get the following message :

  message: Desired error not necessarily achieved due to precision loss.
  success: False
   status: 2
      fun: 26435.40918822449
        x: [ 7.379e+00 -6.967e-02 ...  2.149e-01  4.295e-01]
      nit: 94
      jac: [-4.883e-04  0.000e+00 ... -2.441e-04  4.883e-04]
 hess_inv: [[ 3.566e-04 -7.314e-05 ... -1.261e-04 -1.775e-06]
            [-7.314e-05  5.599e-05 ...  3.216e-05  5.660e-07]
            ...
            [-1.261e-04  3.216e-05 ...  4.981e-04  1.297e-06]
            [-1.775e-06  5.660e-07 ...  1.297e-06  8.380e-06]]
     nfev: 4235
     njev: 121 

The log-likelihood function:

def log_likelihood(params, x, y_lower, y_upper):
    beta = params[:-1]
    sigma = np.abs(params[-1])  

    z_upper = (y_upper - np.dot(x, beta)) / sigma
    z_lower = (y_lower - np.dot(x, beta)) / sigma
    
    epsilon = 1e-10
    ll_interval = np.sum(np.log(np.clip(norm.cdf(z_upper) - norm.cdf(z_lower), epsilon, None)))

    log_likelihood = ll_interval
    
    return -log_likelihood  

and the command that I use to minimize it:

initial_params = np.concatenate([np.zeros(X.shape[1]), [1.0]])
result = minimize(log_likelihood, initial_params, args=(X, y_lower, y_upper), method='BFGS', options={'disp':True})

Does someone has an idea of why the algorithm seems to converge but still giving message that it is not?