Problem 14.1

In [54]:
import math
from __future__ import division

#pade

#solve for b's
def b_coefs(n):
    M = np.zeros((n, n))
    y = np.zeros(n)
    for i in range(n):
        y[i] = -1/math.factorial(n+1+i)
        for j in range(n):
            M[i][j] = 1/math.factorial(n+i-j)
    return np.linalg.pinv(M).dot(y)

#solve for a's
def a_coefs(n, b):
    a = np.zeros(n+1)
    for i in range(n+1):
        a[i] += 1/math.factorial(i)
        for j in range(1, i+1):
            a[i] += b[j-1]/math.factorial(i-j)
    return a

def pade(n,val):
    b = b_coefs(n)
    a = a_coefs(n, b)
    num = a.dot(val**np.arange(0,len(a)))
    denom = 1+b.dot(val**np.arange(1,len(b)+1))
    return num/denom

for n in range(1,6):
    print "["+str(n)+"/"+str(n)+"] pade approximation:"
    val = pade(n, 1)
    print val
    print "error:  " + str(np.abs(val - math.exp(1)))
    print ""

[1/1] pade approximation:
3.0
error:  0.281718171541

[2/2] pade approximation:
2.71428571429
error:  0.00399611417333

[3/3] pade approximation:
2.71830985915
error:  2.80306958853e-05

[4/4] pade approximation:
2.71828171828
error:  1.10177326373e-07

[5/5] pade approximation:
2.71828182874
error:  2.76651146436e-10

Problem 14.2

In [4]:
import numpy as np
import matplotlib.pyplot as plt

x = np.linspace(-10, 10, num=21)
y = np.zeros(21)
for i in range (len(x)):
    if (x[i] > 0): y[i] = 1

        

plt.plot(x, y)
axes = plt.gca()
axes.set_ylim([-0.2,1.2])
plt.show()        

/usr/local/lib/python2.7/site-packages/matplotlib/font_manager.py:273: UserWarning: Matplotlib is building the font cache using fc-list. This may take a moment.
  warnings.warn('Matplotlib is building the font cache using fc-list. This may take a moment.')
In [3]:
#psuedo inverse vandermonde

vandermonde_out_of_sample = [];
x_out_of_sample = np.linspace(-10.5, 10.5, num=22)

for order in [5, 10, 15]:
    n = 21
    A = np.zeros((n, order))

    for i in range(n):
        for j in range(order):
            A[i][j] = x[i]**j

    A_inv = np.linalg.pinv(A)
    coefs = np.matmul(A_inv, y)

    N = 1000
    x_fit = np.linspace(-10, 10, num=N)
    y_fit = np.zeros(N)
    for i in range(N):
        for j in range(order):
            y_fit[i] += coefs[j]*x_fit[i]**j

    out_of_sample = np.zeros(22)
    for i in range(len(x_out_of_sample)):
        for j in range(order):
            out_of_sample[i] += coefs[j]*x_out_of_sample[i]**j
    vandermonde_out_of_sample.append(out_of_sample)

    plt.plot(x_fit, y_fit, label='Order ' + str(order))


plt.plot(x, y)
plt.title('Vandermonde')
axes = plt.gca()
axes.set_ylim([-0.3,1.2])
plt.legend(loc=2)
plt.show()
In [4]:
#RBFs - cubic, uniformly distributed between -10 and 10

rbf_out_of_sample = [];

for order in [5, 10, 15]:
    n = 21
    A = np.zeros((n, order))
    
    centers = np.zeros(order)
    for i in range(order):
        centers[i] = np.random.rand()*20-10

    for i in range(n):
        for j in range(order):
            A[i][j] = np.abs(x[i]-centers[j])**3

    A_inv = np.linalg.pinv(A)
    coefs = np.matmul(A_inv, y)

    N = 1000
    x_fit = np.linspace(-10, 10, num=N)
    y_fit = np.zeros(N)
    for i in range(N):
        for j in range(order):
            y_fit[i] += coefs[j]*np.abs(x_fit[i]-centers[j])**3
            
    out_of_sample = np.zeros(22)
    for i in range(len(x_out_of_sample)):
        for j in range(order):
            out_of_sample[i] += coefs[j]*np.abs(x_out_of_sample[i]-centers[j])**3
    rbf_out_of_sample.append(out_of_sample)

    plt.plot(x_fit, y_fit, label='Order ' + str(order))
    
plt.plot(x, y)
plt.title('RBFs')
axes = plt.gca()
axes.set_ylim([-0.3,1.2])
plt.legend(loc=2)
plt.show() 

#this seems to be highly dependent on the centers of the function
In [21]:
#out of sample error
out_of_sample_gt = np.zeros(22)
for i in range (len(x_out_of_sample)):
    if (x_out_of_sample[i] == 0.5): out_of_sample_gt[i] = 0.5
    elif (x_out_of_sample[i] > 0): out_of_sample_gt[i] = 1

total_error_vdm = []
for i in range(3):
    plt.plot(x_out_of_sample, vandermonde_out_of_sample[i]-out_of_sample_gt)
    total_error_vdm.append(np.sum(np.abs(vandermonde_out_of_sample[i]-out_of_sample_gt)))

print "error:  " + str(total_error_vdm )   

plt.title('out-of-sample error - vandermonde')   
axes = plt.gca()
axes.set_xlim([-10.5,10.5])
plt.show() 

#larger error at higher order

error:  [2.6903613132895727, 2.0511427783907754, 14.214107301343379]
In [23]:
total_error_rbf = []
for i in range(3):
    plt.plot(x_out_of_sample, rbf_out_of_sample[i]-out_of_sample_gt)
    total_error_rbf.append(np.sum(np.abs(rbf_out_of_sample[i]-out_of_sample_gt)))

print "error:  " + str(total_error_rbf )

plt.title('out-of-sample error - rbf')   
axes = plt.gca()
axes.set_xlim([-10.5,10.5])
plt.show() 

#smaller error at higher order

error:  [1.9137117683298182, 1.1594153232054882, 0.55233862941146927]
In [ ]:
order = 10

14.3

In [1]:
x_prev = [0, 0, 0, 1]
x = 0
training_data = []
training_data.append(x_prev)
while (1):
    x = (x_prev[0] + x_prev[3]) % 2
    x_prev.pop()
    training_data.append([x] + x_prev)
    
    x_prev.insert(0,x)
    
    if x_prev == [0, 0, 0, 1]:
        break
training_data.append(x_prev)
        
print training_data

[[0, 0, 0, 1], [1, 0, 0, 0], [1, 1, 0, 0], [1, 1, 1, 0], [1, 1, 1, 1], [0, 1, 1, 1], [1, 0, 1, 1], [0, 1, 0, 1], [1, 0, 1, 0], [1, 1, 0, 1], [0, 1, 1, 0], [0, 0, 1, 1], [1, 0, 0, 1], [0, 1, 0, 0], [0, 0, 1, 0], [0, 0, 0, 1], [0, 0, 0, 1]]
In [11]:
weights = np.ones((4, 4))
weights*=0.25#init with all equal contributions
training_data = np.array(training_data)

c = 0.1

def eval_layer(data, _weights):
    val = _weights.dot(data)
    if (np.abs(val) == val): return 1
    return 0

def eval_entire_layer(data):
    out = []
    for i in range(4):
        out.append(eval_layer(data, weights[i]))
    return np.array(out)

def update_weights(data, output, expected):
    for i in range(4):
        error = expected[i] - output[i]
        weights[i] += error*data*c

#train
for a in range(len(training_data)-1):
    data = training_data[a]
    output = eval_entire_layer(data)
    update_weights(data, output, training_data[a+1])

print weights

[[-0.05 -0.05 -0.05 -0.05]
 [ 0.25 -0.05 -0.05 -0.05]
 [-0.05  0.25 -0.05 -0.05]
 [-0.05 -0.05  0.25  0.05]]
In [12]:
from __future__ import division
percent_correct = 0;
digits_correct = 0

def compare(data, output):
    num = 0
    for i in range(4):
        if (data[i] == output[i]): num+=1
    return num/4

print "no hidden layers"
print ""
print "output       training data"
for a in range(len(training_data)-1):
    data = training_data[a]
    output = eval_entire_layer(data)
    correct = compare(training_data[a+1], output)
    if (correct == 1): percent_correct += 1
    digits_correct += correct
    print str(output)  + "    " + str(training_data[a+1])
    
percent_correct /= len(training_data)-1
digits_correct /= len(training_data)-1
print ""
print "percent correct:  " + str(percent_correct)
print "percent digits correct:  " + str(digits_correct)

no hidden layers

output       training data
[0 0 0 1]    [1 0 0 0]
[0 1 0 0]    [1 1 0 0]
[0 1 1 0]    [1 1 1 0]
[0 1 1 1]    [1 1 1 1]
[0 1 1 1]    [0 1 1 1]
[0 0 1 1]    [1 0 1 1]
[0 1 0 1]    [0 1 0 1]
[0 0 1 0]    [1 0 1 0]
[0 1 0 1]    [1 1 0 1]
[0 1 1 0]    [0 1 1 0]
[0 0 1 1]    [0 0 1 1]
[0 0 0 1]    [1 0 0 1]
[0 1 0 0]    [0 1 0 0]
[0 0 1 0]    [0 0 1 0]
[0 0 0 1]    [0 0 0 1]
[0 0 0 1]    [0 0 0 1]

percent correct:  0.5
percent digits correct:  0.859375