# Using optimization routines from `scipy` and `statsmodels`


```python
%matplotlib inline
```


```python
import scipy.linalg as la
import numpy as np
import scipy.optimize as opt
import matplotlib.pyplot as plt
import pandas as pd
```


```python
np.set_printoptions(precision=3, suppress=True)
```

Using `scipy.optimize`
----

One of the most convenient libraries to use is `scipy.optimize`, since it is already part of the Anaconda installation and it has a fairly intuitive interface.


```python
from scipy import optimize as opt
```

#### Minimizing a univariate function $f: \mathbb{R} \rightarrow \mathbb{R}$


```python
def f(x):
    return x**4 + 3*(x-2)**3 - 15*(x)**2 + 1
```


```python
x = np.linspace(-8, 5, 100)
plt.plot(x, f(x));
```

The [`minimize_scalar`](http://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.minimize_scalar.html#scipy.optimize.minimize_scalar) function will find the minimum, and can also be told to search within given bounds. By default, it uses the Brent algorithm, which combines a bracketing strategy with a parabolic approximation.


```python
opt.minimize_scalar(f, method='Brent')
```




         fun: -803.3955308825884
        nfev: 17
         nit: 11
     success: True
           x: -5.528801125219663




```python
opt.minimize_scalar(f, method='bounded', bounds=[0, 6]) # Checks bounds as well
```




         fun: -54.21003937712762
     message: 'Solution found.'
        nfev: 12
      status: 0
     success: True
           x: 2.668865104039653



### Local and global minima


```python
def f(x, offset):
    return -np.sinc(x-offset)
```


```python
x = np.linspace(-20, 20, 100)
plt.plot(x, f(x, 5));
```


```python
# note how additional function arguments are passed in
sol = opt.minimize_scalar(f, args=(5,))
sol
```




         fun: -0.049029624014074166
        nfev: 15
         nit: 10
     success: True
           x: -1.4843871263953001




```python
plt.plot(x, f(x, 5))
plt.axvline(sol.x, c='red')
pass
```

#### We can try multiple random starts to find the global minimum


```python
lower = np.random.uniform(-20, 20, 100)
upper = lower + 1
sols = [opt.minimize_scalar(f, args=(5,), bracket=(l, u)) for (l, u) in zip(lower, upper)]
```


```python
idx = np.argmin([sol.fun for sol in sols])
sol = sols[idx]
```


```python
plt.plot(x, f(x, 5))
plt.axvline(sol.x, c='red');
```

#### Using a stochastic algorithm

See documentation for the [`basinhopping`](http://docs.scipy.org/doc/scipy-0.14.0/reference/generated/scipy.optimize.basinhopping.html) algorithm, which also works with multivariate scalar optimization. Note that this is heuristic and not guaranteed to find a global minimum.


```python
from scipy.optimize import basinhopping

x0 = 0
sol = basinhopping(f, x0, stepsize=1, minimizer_kwargs={'args': (5,)})
sol
```




                            fun: -1.0
     lowest_optimization_result:       fun: -1.0
     hess_inv: array([[0.30396369]])
          jac: array([1.49011612e-08])
      message: 'Optimization terminated successfully.'
         nfev: 39
          nit: 3
         njev: 13
       status: 0
      success: True
            x: array([5.])
                        message: ['requested number of basinhopping iterations completed successfully']
          minimization_failures: 0
                           nfev: 1668
                            nit: 100
                           njev: 556
                              x: array([5.])




```python
plt.plot(x, f(x, 5))
plt.axvline(sol.x, c='red');
```

### Constrained optimization with `scipy.optimize`

Many real-world optimization problems have constraints - for example, a set of parameters may have to sum to 1.0 (equality constraint), or some parameters may have to be non-negative (inequality constraint). Sometimes, the constraints can be incorporated into the function to be minimized, for example, the non-negativity constraint $p \gt 0$ can be removed by substituting $p = e^q$ and optimizing for $q$. Using such workarounds, it may be possible to convert a constrained optimization problem into an unconstrained one, and use the methods discussed above to solve the problem.

Alternatively, we can use optimization methods that allow the specification of constraints directly in the problem statement as shown in this section. Internally, constraint violation penalties, barriers and Lagrange multipliers are some of the methods used used to handle these constraints. We use the example provided in the Scipy [tutorial](http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html) to illustrate how to set constraints.

We will optimize:

$$
f(x) = -(2xy + 2x - x^2 -2y^2)
$$
subject to the constraint
$$
x^3 - y = 0 \\
y - (x-1)^4 - 2 \ge 0
$$
and the bounds
$$
0.5 \le x \le 1.5 \\
1.5 \le y \le 2.5
$$


```python
def f(x):
    return -(2*x[0]*x[1] + 2*x[0] - x[0]**2 - 2*x[1]**2)
```


```python
x = np.linspace(0, 3, 100)
y = np.linspace(0, 3, 100)
X, Y = np.meshgrid(x, y)
Z = f(np.vstack([X.ravel(), Y.ravel()])).reshape((100,100))
plt.contour(X, Y, Z, np.arange(-1.99,10, 1), cmap='jet');
plt.plot(x, x**3, 'k:', linewidth=1)
plt.plot(x, (x-1)**4+2, 'k:', linewidth=1)
plt.fill([0.5,0.5,1.5,1.5], [2.5,1.5,1.5,2.5], alpha=0.3)
plt.axis([0,3,0,3])
```

To set constraints, we pass in a dictionary with keys `type`, `fun` and `jac`. Note that the inequality constraint assumes a $C_j x \ge 0$ form. As usual, the `jac` is optional and will be numerically estimated if not provided.


```python
cons = ({'type': 'eq',
         'fun' : lambda x: np.array([x[0]**3 - x[1]]),
         'jac' : lambda x: np.array([3.0*(x[0]**2.0), -1.0])},
        {'type': 'ineq',
         'fun' : lambda x: np.array([x[1] - (x[0]-1)**4 - 2])})

bnds = ((0.5, 1.5), (1.5, 2.5))
```


```python
x0 = [0, 2.5]
```

Unconstrained optimization


```python
ux = opt.minimize(f, x0, constraints=None)
ux
```




          fun: -1.9999999999996365
     hess_inv: array([[0.99829115, 0.50105608],
           [0.50105608, 0.4993548 ]])
          jac: array([ 1.25169754e-06, -1.41561031e-06])
      message: 'Optimization terminated successfully.'
         nfev: 24
          nit: 5
         njev: 6
       status: 0
      success: True
            x: array([2.00000043, 0.99999985])



Constrained optimization


```python
cx = opt.minimize(f, x0, bounds=bnds, constraints=cons)
cx
```




         fun: 2.0499154720235886
         jac: array([-3.48747951,  5.49674535])
     message: 'Optimization terminated successfully.'
        nfev: 25
         nit: 6
        njev: 6
      status: 0
     success: True
           x: array([1.26089314, 2.00463288])




```python
x = np.linspace(0, 3, 100)
y = np.linspace(0, 3, 100)
X, Y = np.meshgrid(x, y)
Z = f(np.vstack([X.ravel(), Y.ravel()])).reshape((100,100))
plt.contour(X, Y, Z, np.arange(-1.99,10, 1), cmap='jet');
plt.plot(x, x**3, 'k:', linewidth=1)
plt.plot(x, (x-1)**4+2, 'k:', linewidth=1)
plt.text(ux['x'][0], ux['x'][1], 'x', va='center', ha='center', size=20, color='blue')
plt.text(cx['x'][0], cx['x'][1], 'x', va='center', ha='center', size=20, color='red')
plt.fill([0.5,0.5,1.5,1.5], [2.5,1.5,1.5,2.5], alpha=0.3)
plt.axis([0,3,0,3]);
```

## Some applications of optimization

### Finding paraemeters for ODE models

This is a specialized application of `curve_fit`, in which the curve to be fitted is defined implicitly by an ordinary differential equation 
$$
\frac{dx}{dt} = -kx
$$
and we want to use observed data to estimate the parameters $k$ and the initial value $x_0$. Of course this can be explicitly solved but the same approach can be used to find multiple parameters for $n$-dimensional systems of ODEs.

[A more elaborate example for fitting a system of ODEs to model the zombie apocalypse](http://adventuresinpython.blogspot.com/2012/08/fitting-differential-equation-system-to.html)


```python
from scipy.integrate import odeint

def f(x, t, k):
    """Simple exponential decay."""
    return -k*x

def x(t, k, x0):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """
    x = odeint(f, x0, t, args=(k,))
    return x.ravel()
```


```python
# True parameter values
x0_ = 10
k_ = 0.1*np.pi

# Some random data genererated from closed form solution plus Gaussian noise
ts = np.sort(np.random.uniform(0, 10, 200))
xs = x0_*np.exp(-k_*ts) + np.random.normal(0,0.1,200)

popt, cov = opt.curve_fit(x, ts, xs)
k_opt, x0_opt = popt

print("k = %g" % k_opt)
print("x0 = %g" % x0_opt)
```

    k = 0.31352
    x0 = 9.54446



```python
import matplotlib.pyplot as plt
t = np.linspace(0, 10, 100)
plt.plot(ts, xs, 'r.', t, x(t, k_opt, x0_opt), '-');
```

### Another example of fitting a system of ODEs using the `lmfit` package

You may have to install the [`lmfit`](http://cars9.uchicago.edu/software/python/lmfit/index.html) package using `pip` and restart your kernel. The `lmfit` algorithm is another wrapper around `scipy.optimize.leastsq` but allows for richer model specification and more diagnostics.


```python
#! python3 -m pip install --quiet lmfit
```


```python
from lmfit import minimize, Parameters, Parameter, report_fit
import warnings
```


```python
def f(xs, t, ps):
    """Lotka-Volterra predator-prey model."""
    try:
        a = ps['a'].value
        b = ps['b'].value
        c = ps['c'].value
        d = ps['d'].value
    except:
        a, b, c, d = ps
        
    x, y = xs
    return [a*x - b*x*y, c*x*y - d*y]

def g(t, x0, ps):
    """
    Solution to the ODE x'(t) = f(t,x,k) with initial condition x(0) = x0
    """
    x = odeint(f, x0, t, args=(ps,))
    return x

def residual(ps, ts, data):
    x0 = ps['x0'].value, ps['y0'].value
    model = g(ts, x0, ps)
    return (model - data).ravel()

t = np.linspace(0, 10, 100)
x0 = np.array([1,1])

a, b, c, d = 3,1,1,1
true_params = np.array((a, b, c, d))

np.random.seed(123)
data = g(t, x0, true_params)
data += np.random.normal(size=data.shape)

# set parameters incluing bounds
params = Parameters()
params.add('x0', value= float(data[0, 0]), min=0, max=10)  
params.add('y0', value=float(data[0, 1]), min=0, max=10)  
params.add('a', value=2.0, min=0, max=10)
params.add('b', value=2.0, min=0, max=10)
params.add('c', value=2.0, min=0, max=10)
params.add('d', value=2.0, min=0, max=10)

# fit model and find predicted values
result = minimize(residual, params, args=(t, data), method='leastsq')
final = data + result.residual.reshape(data.shape)

# plot data and fitted curves
plt.plot(t, data, 'o')
plt.plot(t, final, '-', linewidth=2);

# display fitted statistics
report_fit(result)
```

#### Optimization of graph node placement

To show the many different applications of optimization, here is an example using optimization to change the layout of nodes of a graph. We use a physical analogy - nodes are connected by springs, and the springs resist deformation from their natural length $l_{ij}$. Some nodes are pinned to their initial locations while others are free to move. Because the initial configuration of nodes does not have springs at their natural length, there is tension resulting in a high potential energy $U$, given by the physics formula shown below. Optimization finds the configuration of lowest potential energy given that some nodes are fixed (set up as boundary constraints on the positions of the nodes).

$$
U = \frac{1}{2}\sum_{i,j=1}^n ka_{ij}\left(||p_i - p_j||-l_{ij}\right)^2
$$

Note that the ordination algorithm Multi-Dimensional Scaling (MDS) works on a very similar idea - take a high dimensional data set in $\mathbb{R}^n$, and project down to a lower dimension ($\mathbb{R}^k$) such that the sum of distances $d_n(x_i, x_j) - d_k(x_i, x_j)$, where $d_n$ and $d_k$ are some measure of distance between two points $x_i$ and $x_j$ in $n$ and $d$ dimension respectively, is minimized. MDS is often used in exploratory analysis of high-dimensional data to get some intuitive understanding of its "structure".


```python
from scipy.spatial.distance import pdist, squareform
```

- P0 is the initial location of nodes
- P is the minimal energy location of nodes given constraints
- A is a connectivity matrix - there is a spring between $i$ and $j$ if $A_{ij} = 1$
- $L_{ij}$ is the resting length of the spring connecting $i$ and $j$
- In addition, there are a number of `fixed` nodes whose positions are pinned.


```python
n = 20
k = 1 # spring stiffness
P0 = np.random.uniform(0, 5, (n,2)) 
A = np.ones((n, n))
A[np.tril_indices_from(A)] = 0
L = A.copy()
```


```python
L.astype('int')
```




    array([[0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1],
           [0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0]])




```python
def energy(P):
    P = P.reshape((-1, 2))
    D = squareform(pdist(P))
    return 0.5*(k * A * (D - L)**2).sum()
```


```python
D0 = squareform(pdist(P0))
E0 = 0.5* k * A * (D0 - L)**2
```


```python
D0[:5, :5]
```




    array([[0.        , 5.03948582, 0.92078862, 1.75848046, 1.98999961],
           [5.03948582, 0.        , 5.54600931, 3.41449563, 4.9653244 ],
           [0.92078862, 5.54600931, 0.        , 2.13273908, 2.88766898],
           [1.75848046, 3.41449563, 2.13273908, 0.        , 2.76182274],
           [1.98999961, 4.9653244 , 2.88766898, 2.76182274, 0.        ]])




```python
E0[:5, :5]
```




    array([[0.00000000e+00, 8.15872283e+00, 3.13722121e-03, 2.87646305e-01,
            4.90049617e-01],
           [0.00000000e+00, 0.00000000e+00, 1.03331003e+01, 2.91489457e+00,
            7.86189881e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 6.41548915e-01,
            1.78164710e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
            1.55200968e+00],
           [0.00000000e+00, 0.00000000e+00, 0.00000000e+00, 0.00000000e+00,
            0.00000000e+00]])




```python
energy(P0.ravel())
```




    479.13650839857024




```python
# fix the position of the first few nodes just to show constraints
fixed = 4
bounds = (np.repeat(P0[:fixed,:].ravel(), 2).reshape((-1,2)).tolist() + 
          [[None, None]] * (2*(n-fixed)))
bounds[:fixed*2+4]
```




    [[1.191249528562059, 1.191249528562059],
     [4.0389554314507805, 4.0389554314507805],
     [4.474891439430058, 4.474891439430058],
     [0.216114460398234, 0.216114460398234],
     [1.5097341813135952, 1.5097341813135952],
     [4.902910992971438, 4.902910992971438],
     [2.6975241127686767, 2.6975241127686767],
     [3.1315468085492815, 3.1315468085492815],
     [None, None],
     [None, None],
     [None, None],
     [None, None]]




```python
sol = opt.minimize(energy, P0.ravel(), bounds=bounds)
```

#### Visualization

Original placement is BLUE
Optimized arrangement is RED.


```python
plt.scatter(P0[:, 0], P0[:, 1], s=25)
P = sol.x.reshape((-1,2))
plt.scatter(P[:, 0], P[:, 1], edgecolors='red', facecolors='none', s=30, linewidth=2);
```

Optimization of standard statistical models
---

When we solve standard statistical problems, an optimization procedure similar to the ones discussed here is performed. For example, consider multivariate logistic regression - typically, a Newton-like algorithm known as iteratively reweighted least squares (IRLS) is used to find the maximum likelihood estimate for the generalized linear model family. However, using one of the multivariate scalar minimization methods shown above will also work, for example, the BFGS minimization algorithm. 

The take home message is that there is nothing magic going on when Python or R fits a statistical model using a formula - all that is happening is that the objective function is set to be the negative of the log likelihood, and the minimum found using some first or second order optimization algorithm.


```python
import statsmodels.api as sm
```

### Logistic regression as optimization

Suppose we have a binary outcome measure $Y \in {0,1}$ that is conditinal on some input variable (vector) $x \in (-\infty, +\infty)$. Let the conditioanl probability be $p(x) = P(Y=y | X=x)$. Given some data, one simple probability model is $p(x) = \beta_0 + x\cdot\beta$ - i.e. linear regression. This doesn't really work for the obvious reason that $p(x)$ must be between 0 and 1 as $x$ ranges across the real line. One simple way to fix this is to use the transformation $g(x) = \frac{p(x)}{1 - p(x)} = \beta_0 + x.\beta$. Solving for $p$, we get
$$
p(x) = \frac{1}{1 + e^{-(\beta_0 + x\cdot\beta)}}
$$
As you all know very well, this is logistic regression.

Suppose we have $n$ data points $(x_i, y_i)$ where $x_i$ is a vector of features and $y_i$ is an observed class (0 or 1). For each event, we either have "success" ($y = 1$) or "failure" ($Y = 0$), so the likelihood looks like the product of Bernoulli random variables. According to the logistic model, the probability of success is $p(x_i)$ if $y_i = 1$ and $1-p(x_i)$ if $y_i = 0$. So the likelihood is
$$
L(\beta_0, \beta) = \prod_{i=1}^n p(x_i)^y(1-p(x_i))^{1-y}
$$
and the log-likelihood is 
\begin{align}
l(\beta_0, \beta) &= \sum_{i=1}^{n} y_i \log{p(x_i)} + (1-y_i)\log{1-p(x_i)} \\
&= \sum_{i=1}^{n} \log{1-p(x_i)} + \sum_{i=1}^{n} y_i \log{\frac{p(x_i)}{1-p(x_i)}} \\
&= \sum_{i=1}^{n} -\log 1 + e^{\beta_0 + x_i\cdot\beta} + \sum_{i=1}^{n} y_i(\beta_0 + x_i\cdot\beta)
\end{align}

Using the standard 'trick', if we augment the matrix $X$ with a column of 1s, we can write $\beta_0 + x_i\cdot\beta$ as just $X\beta$.


```python
df_ = pd.read_csv("binary.csv")
df_.columns = df_.columns.str.lower()
df_.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>admit</th>
      <th>gre</th>
      <th>gpa</th>
      <th>rank</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>380</td>
      <td>3.61</td>
      <td>3</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>660</td>
      <td>3.67</td>
      <td>3</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>800</td>
      <td>4.00</td>
      <td>1</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>640</td>
      <td>3.19</td>
      <td>4</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0</td>
      <td>520</td>
      <td>2.93</td>
      <td>4</td>
    </tr>
  </tbody>
</table>
</div>




```python
# We will ignore the rank categorical value

cols_to_keep = ['admit', 'gre', 'gpa']
df = df_[cols_to_keep]
df.insert(1, 'dummy', 1)
df.head()
```




<div>
<style scoped>
    .dataframe tbody tr th:only-of-type {
        vertical-align: middle;
    }

    .dataframe tbody tr th {
        vertical-align: top;
    }

    .dataframe thead th {
        text-align: right;
    }
</style>
<table border="1" class="dataframe">
  <thead>
    <tr style="text-align: right;">
      <th></th>
      <th>admit</th>
      <th>dummy</th>
      <th>gre</th>
      <th>gpa</th>
    </tr>
  </thead>
  <tbody>
    <tr>
      <th>0</th>
      <td>0</td>
      <td>1</td>
      <td>380</td>
      <td>3.61</td>
    </tr>
    <tr>
      <th>1</th>
      <td>1</td>
      <td>1</td>
      <td>660</td>
      <td>3.67</td>
    </tr>
    <tr>
      <th>2</th>
      <td>1</td>
      <td>1</td>
      <td>800</td>
      <td>4.00</td>
    </tr>
    <tr>
      <th>3</th>
      <td>1</td>
      <td>1</td>
      <td>640</td>
      <td>3.19</td>
    </tr>
    <tr>
      <th>4</th>
      <td>0</td>
      <td>1</td>
      <td>520</td>
      <td>2.93</td>
    </tr>
  </tbody>
</table>
</div>



### Solving as a GLM with IRLS

This is very similar to what you would do in R, only using Python's `statsmodels` package. The GLM solver uses a special variant of Newton's method known as iteratively reweighted least squares (IRLS), which will be further desribed in the lecture on multivarite and constrained optimizaiton.


```python
model = sm.GLM.from_formula('admit ~ gre + gpa', 
                            data=df, family=sm.families.Binomial())
fit = model.fit()
fit.summary()
```




<table class="simpletable">
<caption>Generalized Linear Model Regression Results</caption>
<tr>
  <th>Dep. Variable:</th>         <td>admit</td>      <th>  No. Observations:  </th>  <td>   400</td> 
</tr>
<tr>
  <th>Model:</th>                  <td>GLM</td>       <th>  Df Residuals:      </th>  <td>   397</td> 
</tr>
<tr>
  <th>Model Family:</th>        <td>Binomial</td>     <th>  Df Model:          </th>  <td>     2</td> 
</tr>
<tr>
  <th>Link Function:</th>         <td>logit</td>      <th>  Scale:             </th> <td>  1.0000</td>
</tr>
<tr>
  <th>Method:</th>                <td>IRLS</td>       <th>  Log-Likelihood:    </th> <td> -240.17</td>
</tr>
<tr>
  <th>Date:</th>            <td>Tue, 06 Apr 2021</td> <th>  Deviance:          </th> <td>  480.34</td>
</tr>
<tr>
  <th>Time:</th>                <td>22:23:02</td>     <th>  Pearson chi2:      </th>  <td>  398.</td> 
</tr>
<tr>
  <th>No. Iterations:</th>          <td>4</td>        <th>                     </th>     <td> </td>   
</tr>
<tr>
  <th>Covariance Type:</th>     <td>nonrobust</td>    <th>                     </th>     <td> </td>   
</tr>
</table>
<table class="simpletable">
<tr>
      <td></td>         <th>coef</th>     <th>std err</th>      <th>z</th>      <th>P>|z|</th>  <th>[0.025</th>    <th>0.975]</th>  
</tr>
<tr>
  <th>Intercept</th> <td>   -4.9494</td> <td>    1.075</td> <td>   -4.604</td> <td> 0.000</td> <td>   -7.057</td> <td>   -2.842</td>
</tr>
<tr>
  <th>gre</th>       <td>    0.0027</td> <td>    0.001</td> <td>    2.544</td> <td> 0.011</td> <td>    0.001</td> <td>    0.005</td>
</tr>
<tr>
  <th>gpa</th>       <td>    0.7547</td> <td>    0.320</td> <td>    2.361</td> <td> 0.018</td> <td>    0.128</td> <td>    1.381</td>
</tr>
</table>



### Or use R


```python
%load_ext rpy2.ipython
```


```r
%%R -i df
m <- glm(admit ~ gre + gpa, data=df, family="binomial")
summary(m)
```

    /opt/conda/lib/python3.6/site-packages/rpy2/robjects/pandas2ri.py:191: FutureWarning: from_items is deprecated. Please use DataFrame.from_dict(dict(items), ...) instead. DataFrame.from_dict(OrderedDict(items)) may be used to preserve the key order.
      res = PandasDataFrame.from_items(items)



    
    Call:
    glm(formula = admit ~ gre + gpa, family = "binomial", data = df)
    
    Deviance Residuals: 
        Min       1Q   Median       3Q      Max  
    -1.2730  -0.8988  -0.7206   1.3013   2.0620  
    
    Coefficients:
                 Estimate Std. Error z value Pr(>|z|)    
    (Intercept) -4.949378   1.075093  -4.604 4.15e-06 ***
    gre          0.002691   0.001057   2.544   0.0109 *  
    gpa          0.754687   0.319586   2.361   0.0182 *  
    ---
    Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
    
    (Dispersion parameter for binomial family taken to be 1)
    
        Null deviance: 499.98  on 399  degrees of freedom
    Residual deviance: 480.34  on 397  degrees of freedom
    AIC: 486.34
    
    Number of Fisher Scoring iterations: 4
    



### Home-brew logistic regression using a generic minimization function

This is to show that there is no magic going on - you can write the function to minimize directly from the log-likelihood equation and run a minimizer. It will be more accurate if you also provide the derivative (+/- the Hessian for second order methods), but using just the function and numerical approximations to the derivative will also work. As usual, this is for illustration so you understand what is going on - when there is a library function available, you should probably use that instead.


```python
def f(beta, y, x):
    """Minus log likelihood function for logistic regression."""
    return -((-np.log(1 + np.exp(np.dot(x, beta)))).sum() + (y*(np.dot(x, beta))).sum())
```


```python
beta0 = np.zeros(3)
opt.minimize(f, beta0, args=(df['admit'], df.loc[:, 'dummy':]), method='BFGS', options={'gtol':1e-2})
```




          fun: 240.17199014098995
     hess_inv: array([[ 1.12219014e+00, -2.92706683e-04, -2.70998572e-01],
           [-2.92706683e-04,  1.11080617e-06, -1.10402280e-04],
           [-2.70998572e-01, -1.10402280e-04,  9.81365109e-02]])
          jac: array([-3.81469727e-06, -2.06756592e-03, -9.53674316e-06])
      message: 'Optimization terminated successfully.'
         nfev: 65
          nit: 8
         njev: 13
       status: 0
      success: True
            x: array([-4.94931338e+00,  2.69042448e-03,  7.54713547e-01])



### Optimization with `sklearn`

There are also many optimization routines in the `scikit-learn` package, as you already know from the previous lectures. Many machine learning problems essentially boil down to the minimization of some appropriate loss function.

### Resources

- [Scipy Optimize reference](http://docs.scipy.org/doc/scipy/reference/optimize.html)
- [Scipy Optimize tutorial](http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html)
- [LMFit - a modeling interface for nonlinear least squares problems](http://cars9.uchicago.edu/software/python/lmfit/index.html)
- [CVXpy- a modeling interface for convex optimization problems](https://github.com/cvxgrp/cvxpy)
- [Quasi-Newton methods](http://en.wikipedia.org/wiki/Quasi-Newton_method)
- [Convex optimization book by Boyd & Vandenberghe](http://stanford.edu/~boyd/cvxbook/)
- [Nocedal and Wright textbook](http://www.springer.com/us/book/9780387303031)

66cdf962e16fb866cc61741c60030b597fc4cf13
ipynb
Jupyter Notebook
notebooks/copies/lectures/T07D_Optimization_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/T07D_Optimization_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
notebooks/copies/lectures/T07D_Optimization_Examples.ipynb
robkravec/sta-663-2021
4dc8018f7b172eaf81da9edc33174768ff157939
MIT
Qwen/Qwen-72B
1. YES
2. YES