<jupyter_start><jupyter_text>Optimization<jupyter_code>!pip install numdifftools
import numpy as np
import timeit
from scipy import optimize as opt
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d.axes3d import get_test_data
import numdifftools as nd<jupyter_output><empty_output><jupyter_text>Scalar functions<jupyter_code>def f(x):
    return (x+3)*(x-1)**2
X = np.arange(-5, 5, .1)
Y = f(X)
plt.plot(X, Y)
plt.show()
%timeit opt.minimize_scalar(f, bracket=(-2, 2, 4), method='golden')
%timeit opt.minimize_scalar(f, bracket=(-2, 2, 4), method='brent')<jupyter_output>714 µs ± 33.8 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)<jupyter_text>Function to calculate first derivative of $f$ at $x_{0}$<jupyter_code>def d1(x, f, h=1e-5):
    return (f(x+h)-f(x-h))/(2*h)
# Test
d1(5, f)<jupyter_output><empty_output><jupyter_text>Function to calculate second derivative of $f$ at $x_{0}$<jupyter_code>def d2(x, f, h=1e-5):
    return (f(x+2*h)+f(x-2*h)-2*f(x))/(2*h)**2
# Test
d2(5, f)<jupyter_output><empty_output><jupyter_text>Newton-Raphson<jupyter_code>def my_newton(f, x0, h=1e-5, delta=1e-5, eps=1e-5):
    
    dif = 1
    der = 1
    i = 0
    
    while dif>eps or der>delta:
        
        i += 1
        d = d1(x0, f)
        x1 = x0 - (d/d2(x0, f))
        dif = abs(x1-x0)
        der = abs(d)
        x0 = x1
        if i > 100:
            break
    
    return x1
%timeit my_newton(f, 10)<jupyter_output>74.6 µs ± 4.27 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)<jupyter_text>Derivatives<jupyter_code>df1 = nd.Derivative(f)
df2 = nd.Derivative(f, n=2)
%timeit opt.minimize(f, 10, method='Newton-CG', jac=df1, hess=df2)<jupyter_output>185 ms ± 10.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)<jupyter_text>Refresher of matrix algebra Quadratic form: $Q(x_{1}, x_{2}) = a_{1}\,x_{1}^2 + a_{2}\,x_{2}^2 + a_{3}\,x_{1}\,x_{2}$ This a algebraic expression can be written with matrix notation:$   x=  \left[ {\begin{array}{c}   x_{1} \\   x_{2} \\  \end{array} } \right]$ and  $   A=  \left[ {\begin{array}{c}   a_{1} & a_{3}/2 \\   a_{3}/2 & a_{2} \\  \end{array} } \right]$ imply $Q(x)=x^{T}\,A\,x$. Notice that $A$ is a square, **symmetric** matrix. Therefore, there is a one-to-one relationship between quadratic forms and symmetric matrices. *Definitions*: * A quadratic form is *positive definite* iff $Q(x)>0$ for any $x\neq0$.* A quadratic form is *negative definite* iff $Q(x)<0$ for any $x\neq0$.<jupyter_code>%matplotlib notebook
plt.rcParams['figure.figsize'] = [8, 3]
fig = plt.figure()
fig.suptitle('Definite Quadratic Forms', fontsize=16)

x = np.linspace(-5, 5, 101)
y = np.linspace(-5, 5, 101)
X, Y = np.meshgrid(x, y)

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_wireframe(X, Y, X**2+Y**2, rstride=10, cstride=10)
ax1.set_title('Positive Definite')

ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_wireframe(X, Y, -X**2-Y**2, rstride=10, cstride=10)
ax2.set_title('Negative Definite')

plt.show()<jupyter_output><empty_output><jupyter_text>*Definitions*: * A quadratic form is *positive semi-definite* iff $Q(x)\geq0$ for any $x$.* A quadratic form is *negative semi-definite* iff $Q(x)\leq0$ for any $x$.<jupyter_code>%matplotlib notebook
fig = plt.figure()
fig.suptitle('Semi-Definite Quadratic Forms', fontsize=16)

ax1 = fig.add_subplot(121, projection='3d')
ax1.plot_wireframe(X, Y, (X+Y)**2, rstride=10, cstride=10)
ax1.set_title('Positive Semi-Definite')

ax2 = fig.add_subplot(122, projection='3d')
ax2.plot_wireframe(X, Y, -(X+Y)**2, rstride=10, cstride=10)
ax2.set_title('Negative Semi-Definite')

plt.show()<jupyter_output><empty_output><jupyter_text>*Definitions*: * A quadratic form is *indefinite* iff $Q(x)>0$ for some $x$ and $Q(x)<0$ for another $x$.<jupyter_code>%matplotlib notebook
plt.rcParams['figure.figsize'] = [5, 3]
fig = plt.figure()
fig.suptitle('Indefinite Quadratic Forms', fontsize=16)

ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, X**2-Y**2, rstride=10, cstride=10)

plt.show()<jupyter_output><empty_output><jupyter_text>Testing definiteness<jupyter_code>def principal_minors(M):
    if M.ndim == 2 and np.all(M == M.transpose()):
        lpm = []
        for k in np.arange(0, M.shape[0], 1):
            lpm.append(np.linalg.det(M[0:k+1, 0:k+1]))
        return lpm
A = np.array([[1, 2, 3], [2, 3, 4], [3, 4, 9]])
principal_minors(A)
def positive_definite(M):
    if M.ndim == 2 and np.all(M == M.transpose()):
        return np.all(np.array(principal_minors(M))>0)
positive_definite(A)
def negative_definite(M):
    if M.ndim == 2 and np.all(M == M.transpose()):
        n = M.shape[0]
        c = [(-1)**(i+1) for i in np.arange(0, n, 1)]
        return np.all((np.array(c)*np.array(principal_minors(M)))>0)
negative_definite(A)
def indefinite(M):
    if M.ndim == 2 and np.all(M == M.transpose()):
        if not positive_definite(M) and not negative_definite(M):
            n = M.shape[0]
            return np.all(np.array(principal_minors(M)) != 0)
indefinite(A)<jupyter_output><empty_output><jupyter_text>Eigenvalues, Eigenvectors and Eigen Decomposition$A$: square matrix ($n \times n$), $x$: column vector ($n \times 1$), $\lambda$: scalar. *Definition*: $\lambda$ is an **eigenvalue** of $A$ with corresponding $x$ **eigenvector** iff $A\,x=\lambda\,x$.Eigenvalues must solve $\mathrm{det}(A-\lambda\,I)=0$.*Result*: If $A$ is a real symmetric matrix, all of its eigenvalues are real numbers, and eigenvectors corresponding to distinct eigenvalues are orthogonal.*Eigen Decomposition of Symmetric Matrices*: Let $P = [x_1, x_2, ..., x_n]$ and $D=diag(\lambda_1, \lambda_2, ..., \lambda_n)$. Then $A = P\,D\,P^{-1}$.If $\lambda_i>0, \forall i$, then $A$ is **positive definite**. If $\lambda_iIf $\lambda_i\geq0, \forall i$, then $A$ is **positive semi-definite**. If $\lambda_i\leq0, \forall i$, then $A$ is **negative semi-definite**.<jupyter_code>np.linalg.eig(A)
def pos_def(A):
    return (np.all(np.linalg.eig(A)[0]>0))
pos_def(A)<jupyter_output><empty_output><jupyter_text>Multivariate functions<jupyter_code>def f(x):
    return x[0]**2+x[1]**2-50*np.sin(.5*(x[0]+x[1]))
x = np.array([1, 2, 3])
y = np.array([20, 40, 60])
X, Y = np.meshgrid(x, y)
X
Y
%matplotlib notebook
plt.rcParams['figure.figsize'] = [5, 3]

x = np.linspace(-10, 10, 101)
y = np.linspace(-10, 10, 101)
X, Y = np.meshgrid(x, y)

fig = plt.figure()
fig.suptitle('$x^2+y^2-50\,sin(.5(x+y))$', fontsize=16)

ax = fig.add_subplot(111, projection='3d')
ax.plot_wireframe(X, Y, f([X, Y]), rstride=10, cstride=10)

plt.show()
%matplotlib notebook
plt.rcParams['figure.figsize'] = [5, 3]
x = np.linspace(-10, 10, 101)
y = np.linspace(-10, 10, 101)
X, Y = np.meshgrid(x, y)

fig = plt.figure()

ax = fig.add_subplot(111)
ax.contour(X, Y, f([X, Y]))

plt.show()
import numpy as np
from scipy.optimize import minimize<jupyter_output><empty_output><jupyter_text>Nelder-Mead<jupyter_code>x0 = np.array([0, 0])
res = minimize(f, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True})
print(res.x)
f(res.x)
res<jupyter_output><empty_output><jupyter_text>Gradient and Hessian<jupyter_code>fgrad = nd.Gradient(f)
fhess = nd.Hessian(f)
type(fgrad)<jupyter_output><empty_output><jupyter_text>Conjugate Gradient<jupyter_code>x0 = np.array([0, 0])
res = minimize(f, x0, method='CG', jac=fgrad, options={'disp': True})
print(res.x)<jupyter_output>Optimization terminated successfully.
         Current function value: -45.431122
         Iterations: 3
         Function evaluations: 8
         Gradient evaluations: 8
[1.45419657 1.45419657]<jupyter_text>Broyden, Fletcher, Goldfarb, and Shanno (BFGS)<jupyter_code>x0 = np.array([0, 0])
res = minimize(f, x0, method='BFGS', jac=fgrad, options={'gtol': 1e-8, 'disp': True})
print(res.x)
f(res.x)<jupyter_output>Optimization terminated successfully.
         Current function value: -45.431122
         Iterations: 5
         Function evaluations: 7
         Gradient evaluations: 7
[1.45419657 1.45419657]<jupyter_text>Newton-CG<jupyter_code>x0 = np.array([0, 0])
res = minimize(f, x0, method='Newton-CG', jac=fgrad, hess=fhess, options={'xtol': 1e-8, 'disp': True})
print(res.x)
f(res.x)<jupyter_output>Optimization terminated successfully.
         Current function value: -45.431122
         Iterations: 5
         Function evaluations: 9
         Gradient evaluations: 13
         Hessian evaluations: 5
[1.45419658 1.45419658]<jupyter_text>Maximum Likelihood<jupyter_code>mu = 18
sigma = 5
n = 30
S = np.random.normal(loc=mu, scale=sigma, size=n)
def likelihood(par):
    # par[0]=mu
    # par[1]=sigma
    
    v = 0
    
    for x in S:
        v += -.5*((x-par[0])/par[1])**2-np.log(par[1])
    
    return -v
x0 = np.array([1, 1])
res = minimize(likelihood, x0, method='nelder-mead', options={'xtol': 1e-8, 'disp': True})
res.x
S.mean()<jupyter_output><empty_output><jupyter_text>Consumer´s problem with CES utility<jupyter_code>def utility(x, *args):
    return -(x**.25+(args[0]-args[1]*x)**.25)
def demand(m, p):
    x = opt.minimize(utility, .3, args=(m, p), method="CG").x[0]
    return (x, m-p*x)
demand(5, 2)
price = np.arange(.1, 5, .1)<jupyter_output><empty_output><jupyter_text>Regressão Múltipla<jupyter_code>x1 = np.random.uniform(0, 100, size=20000)
x2 = np.random.uniform(0, 100, size=20000)
x2 = .9*x2+.1*x1
sigma = 5
y = 10+.4*x1+.7*x2+np.random.normal(loc=0, scale=sigma, size=20000)
def my_ols(beta):  
    return np.sum((y-beta[0]-beta[1]*x1-beta[2]*x2)**2)   
opt.minimize(my_ols, (10.1, .5, .6), method="CG")<jupyter_output><empty_output><jupyter_text>OLS with matrices<jupyter_code>import io
DATA = np.loadtxt("data/regression.csv", delimiter=";")
DATA
DATA.shape
y = DATA[:, 2].reshape(500, 1)
X = DATA[:, 0:2].reshape(500, 2)
c = np.ones(500).reshape(500, 1)
X = np.concatenate((c, X), axis=1)
def my_ols(y, X):
    
    results = {}
    
    (n, k) = X.shape
    results['n'] = n
    results['k'] = k
    
    beta = np.linalg.inv(X.T @ X) @ X.T @ y
    results['beta'] = beta
    
    e = y - X @ beta
    results['residuas'] = e
    
    SSE = sum([e**2 for e in e])[0]
    results['SSE'] = SSE
    
    sigma = np.sqrt(SSE/(n-k))
    results['sigma'] = sigma
    
    V = sigma**2*np.linalg.inv(X.T @ X)
    results['V'] = V
    
    stderr_beta = np.sqrt(np.diag(V))
    results['stderr_beta'] = stderr_beta
    
    t = beta/stderr_beta
    results['t'] = t
    
    # results['p-values'] = 
    
    ybar = np.mean(y)
    SQT = sum([(q-ybar)**2 for q in y])[0]
    results['SQT'] = SQT
    
    results['R2'] = 1-(SSE/SQT)
    
    return results 
res = my_ols(y, X)
res
res['R2']<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18347932560987096, 'fraction_numerical': 0.05010831685033437, 'mean_word_length': 4.161886242100146, 'pattern_counts': {'":': 0, '<': 81, '<?xml version=': 0, '>': 88, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '21', 'dir': 'jupyter-structured-clean-dedup', 'id': '207353', 'license': 'MIT', 'n_tokens_mistral': 4767, 'n_tokens_neox': 4380, 'n_words': 951, 'path': '00-pre-requisitos/2-math/otimização-I.ipynb', 'repo_name': 'sn3fru/datascience_course', 'size': '10617'}