
import numpy as np
import matplotlib.pyplot as plt


def f(x):
    return np.sin(2.*np.pi*x) + x**2 + 2.*(np.exp(-30.*(x-0.6)**2) - np.exp(-30.*((x-0.6)/3.)**2))


xx = np.linspace(0, 1, 1000)
plt.plot(xx, f(xx));


b = np.arange(0, 1, 0.1)
print(b)

[0.  0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9]


v = f(b)
print(v)

[-0.60234762 -0.270305   -0.16577643 -0.3061689  -0.40017296 -0.20279576
 -0.22778525 -0.91385228 -1.45901473 -1.12501067]


# A piecewise constant function
def _G(x, v, b):
    g = v[0]
    idx = 1
    while idx<len(b) and x>b[idx]:
        #print(idx)
        g = v[idx]
        idx +=1
    return g

def G(x, v, b):
    '''
     y = G(x, v, b)
     
     Samples the piecewise constant function defined by the
     values in v at breakpoints in b.
     
     Inputs:
      x   array of x-values
      v   array of constant piece values
      b   array of breakpoints (the left of each piece)
          Note: len(v)=len(b)
    '''
    val = [_G(y, v, b) for y in x]
    return np.array(val)


b = np.linspace(0,1,10)
plt.plot(xx, f(xx))
plt.plot(xx, G(xx, f(b), b), 'r');
plt.plot(b, f(b), 'ko')
max_diff = np.max(abs(G(xx, f(b), b)-f(xx)))
plt.title('Max error = '+str(max_diff));


b = np.linspace(0,1,40)
plt.plot(xx, f(xx))
plt.plot(xx, G(xx, f(b), b), 'r');
max_diff = np.max(abs(G(xx, f(b), b)-f(xx)))
plt.title('Max error = '+str(max_diff));


def find_interval(f, x, eps):
    '''
     delta = find_interval(f, x, eps)
     Find the step delta s.t.
      | f(x) - f(x+delta) | < eps
    '''
    delta = 1.
    v = f(x)
    err = np.max(np.abs(v-f(np.linspace(x, x+delta, 100))))
    while err >= eps:
        delta *= 0.98
        err = np.max(np.abs(v-f(np.linspace(x, x+delta, 100))))
    
    return delta

def build_approx(f, domain, eps):
    '''
     v, b = build_approx(f, domain, eps)
     Find a piecewise constant approximation to f that is
     within eps of f at all points.
    '''
    x = domain[0]  # left side
    b = [x]
    v = [f(x)]
    while x<domain[1]:
        b.append(b[-1]+find_interval(f, x, eps))
        v.append(f(b[-1]))
        x = b[-1]
    b[-1] = domain[1]
    v[-1] = f(b[-1])
    return v, b


eps = 0.1
v, b = build_approx(f, [-1, 1], eps)
print(f'{len(b)} pieces')

72 pieces


plt.plot(xx, f(xx))
plt.plot(xx, G(xx, v, b), 'r');
plt.fill_between(xx, G(xx, v, b)-eps, G(xx, v, b)+eps, color=(1,0,0,0.1))
#plt.plot(xx, G(xx, v, b)+eps, color=(1,0,0,0.2));
#plt.plot(xx, G(xx, v, b)-eps, color=(1,0,0,0.2));
max_diff = np.max(abs(G(xx, v, b)-f(xx)))
plt.title('Max error = '+str(max_diff));


xx = np.linspace(-1,1,5000)
inv_cubic = (lambda x: np.sign(x)*abs(x)**(1./3))
plt.plot(xx, inv_cubic(xx), 'r');


eps = 0.1
v, b = build_approx(inv_cubic, [xx[0], xx[-1]], eps)
print(f'{len(b)} pieces')

22 pieces


plt.plot(xx, inv_cubic(xx))
plt.plot(xx, G(xx, v, b), 'r');
plt.fill_between(xx, G(xx, v, b)-eps, G(xx, v, b)+eps, color=(1,0,0,0.1))
max_diff = np.max(abs(G(xx, v, b)-inv_cubic(xx)))
plt.title('Max error = '+str(max_diff));


plt.hist(abs(np.diff(b)), bins=10);
plt.xlabel(r'$\Delta x$'); plt.ylabel('Frequency');


def logistic(x, w, b):
    g = 1./(1.+np.exp(w*(b-np.array(x))))
    return g


xx = np.linspace(0, 1, 400)
plt.plot(xx, logistic(xx, 10, 0.4));


N = 10
A = np.zeros([len(xx), N])
for idx in range(N):
    b = np.random.random()
    w = np.random.random()*20 - 10.
    A[:, idx] = logistic(xx, w, b)


plt.plot(xx, A);


plt.plot(xx, f(xx));


c = np.linalg.lstsq(A, f(xx), rcond=-1)[0]


f_est = A @ c
f_diag = A @ np.diag(c)
plt.plot(xx, f(xx), '--', lw=3);
plt.plot(xx, f_est,);
plt.legend(['Original', 'Logistic Fit']);

c

array([-7.88271672e+02,  7.41624011e+03,  5.59652611e+01, -3.90384194e+02,
        6.49213904e+02, -1.83925862e+02,  5.17678449e+04, -5.19717314e+04,
        8.24069068e+00, -1.87412486e+02])

Universal Approximation Theorem¶

Complicated Function to Model¶

Split domain into Intervals¶

Fit $f(x)$ with a piecewise constant function¶

Adaptive partitioning¶

A trickier function¶

Logistic Regression¶