title: Estimating logistic regression using BFGS optimization algorithm
body: <p>I tried to use <code>scipy.optimize.minimum</code> to estimate parameters in logistic regression. Before this, I wrote log likelihood function and gradient of log likelihood function. I then used <strong>Nelder-Mead</strong> and <strong>BFGS</strong> algorithm, respectively. <em>Turned out the latter one failed but the former one succeeded.</em> Because <strong>BFGS</strong> algorithm also uses the gradient, I doubt the gradient part (<code>function log_likelihood_gradient(x, y)</code>) suffers from bugs. But it's also likely due to any numerical issue associated with <strong>Nelder-Mead</strong> since I got the warning <code>Desired error not necessarily achieved due to precision loss.</code></p>
  <hr />
  <p>I'm using UCLA's tutorial and dataset (see <a href="http://www.ats.ucla.edu/stat/r/dae/logit.htm" rel="nofollow noreferrer">link</a>) and the correct estimates are,</p>
  <p><a href="https://i.stack.imgur.com/4BnP2m.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/4BnP2m.png" alt="enter image description here" /></a></p>
  <p>FYI, I also highly suggest you reading this <a href="https://stats.stackexchange.com/questions/203740/logistic-regression-scikit-learn-vs-statsmodels">post</a> to understand the problem setup and preprocess on the dataset.</p>
  <hr />
  <p>Now is my code, the first part is data preparation. You can simply run it to get yourself a copy,</p>
  <pre><code>import numpy as np
  from scipy.optimize import minimize
  from patsy import dmatrices
  
  from urllib.request import urlretrieve
  url = 'http://www.ats.ucla.edu/stat/data/binary.csv'
  urlretrieve(url, './admit_rate.csv')
  data = pd.read_csv('./admit_rate.csv')
  
  y, x = dmatrices('admit ~ gre + gpa + C(rank)', data, return_type = 'dataframe')
  y = np.array(y, dtype=np.float64)
  x = np.array(x, dtype=np.float64)
  </code></pre>
  <hr />
  <p>The second part is the MLE part, <code>log_likelihood</code> defines the sum of log-likelihood function, and <code>log_likelihood_gradient</code> defines the gradient of the sum of log-likelihood function. Both come from the formula as below,
  <a href="https://i.stack.imgur.com/APPCp.png" rel="nofollow noreferrer"><img src="https://i.stack.imgur.com/APPCp.png" alt="enter image description here" /></a></p>
  <pre><code>def sigmoid(x):
      &quot;&quot;&quot;
      Logistic function
      &quot;&quot;&quot;
      return 1.0 / (1 + np.exp(-x))
  
  def sigmoid_derivative(x):
      &quot;&quot;&quot;
      Logistic function's first derivative
      &quot;&quot;&quot;
      return sigmoid(x) * (1 - sigmoid(x))
  
  def log_likelihood(x, y):
      &quot;&quot;&quot;
      The sum of log likelihood functions
      &quot;&quot;&quot;
      n1, n2 = x.shape
      def log_likelihood_p(p):
          &quot;&quot;&quot;
          p: betas/parameters in logistic regression, to be estimated
          &quot;&quot;&quot;
          p = p.reshape(n2, 1)
          sig = sigmoid(np.dot(x, p)).reshape(n1, 1)
          return np.sum(y * np.log(sig) + (1 - y) * np.log(1 - sig)) 
      return log_likelihood_p
   
  def log_likelihood_gradient(x, y):
      &quot;&quot;&quot;
      The gradient of the sum of log likelihood functions used in optimization
      &quot;&quot;&quot;
      n1, n2 = x.shape
      def log_likelihood_gradient_p(p):
          &quot;&quot;&quot;
          p: betas/parameters in logistic regression, to be estimated
          &quot;&quot;&quot;
          p = p.reshape(n2, 1)
          xp = np.dot(x, p)
          sig = sigmoid(xp)
          sig_der = sigmoid_derivative(xp)
          return np.sum(y * sig_der / sig * x - (1 - y) * sig_der / (1 - sig) * x, axis=0)
      return log_likelihood_gradient_p
  
  def negate(f):  
      return lambda *args, **kwargs: -f(*args, **kwargs)
  </code></pre>
  <hr />
  <p>The third part is the optimization part, I started with inital value</p>
  <p><em><strong>x0 = np.array([-1, -0.5, -1, -1.5, 0, 1])</strong></em></p>
  <p>for <strong>Nelder-Mead</strong>, it works and outputs</p>
  <p><em><strong>array([ -3.99005691e+00,  -6.75402757e-01,  -1.34021420e+00, -1.55145450e+00,   2.26447746e-03,   8.04046227e-01]</strong></em>.</p>
  <p>But then when I tried <strong>Nelder-Mead</strong> algorithm, it failed at</p>
  <p><em><strong>x0 = np.array([-1, -0.5, -1, -1.5, 0, 1])</strong></em>.</p>
  <p>Even when I gave even closer initial value like the output from <strong>Nelder-Mead</strong> algorithm, it failed with warning <code>Desired error not necessarily achieved due to precision loss</code>.</p>
  <p><strong>Nelder-Mead</strong></p>
  <pre><code>x0 = np.array([-1, -0.5, -1, -1.5, 0, 1])
  estimator1 = minimize(negate(log_likelihood(x, y)), x0, 
                 method='nelder-mead', options={'disp': True})
  print(estimator1.success)
  estimator1.x
  </code></pre>
  <p><strong>BFGS</strong></p>
  <pre><code>x0 = np.array([ -3.99005691e+00,  -6.75402757e-01,  -1.34021420e+00,
          -1.55145450e+00,   2.26447746e-03,   8.04046227e-01])
  estimator2 = minimize(negate(log_likelihood(x, y)), x0, method='BFGS', 
                       jac=log_likelihood_gradient(x, y), options={'disp': True})
  print(estimator2.success)
  </code></pre>
  <hr />
  <p>I know my post is a bit long, but I'd glad to answer for any questions you throw to me and can anyone help?</p>