Previously I tried

opt1 = minimize(logit, args = (df), x0=x_start,method = 'L-BFGS-B')
hessinv1 = opt1.hess_inv.todense()
hessinvm1 = np.asmatrix(hessinv1)
diag1 = np.sqrt(np.diagonal(hessinvm1))

However, it turns out that this Hessian matrix is very unreliable and the results of standard errors are not correct, which has also been mentioned here (LBFGS: Accuracy of Hessian approximation). Is there any way that we could obtain the standard error after the optimization process?