Problem 12.1

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

N = 100
stdev = 0.5
x = np.random.rand(N)
y = 2+3*x
for i in range(N):
    y[i] += np.random.normal(0, stdev)

actual,= plt.plot(np.linspace(0,1,N), 2+3*np.linspace(0,1,N))
plt.scatter(x, y)
axes = plt.gca()
axes.set_xlim([0,1])
plt.legend([actual], ['y = 2 + 3x'], loc=2)
plt.show()

SVD

In [201]:
def SVD(_x, _y):
    A = []
    A = np.matrix(np.append(np.ones([N,1]), _x.reshape(N,1), 1))
    U, s, V = np.linalg.svd(A)

    #invert s
    s_inv = np.copy(s);
    for i in range(len(s_inv)):
        if (s_inv[i] == 0):
            s_inv[i] = 0;
        else :
            s_inv[i] = 1/s_inv[i];

    w = np.zeros((2,N))
    w[:2,:2] = np.diag(s_inv)

    # print (V.T).shape
    # print w.shape
    # print (U.T).shape
    # print (y.T).shape

    c = V.T*w*U.T*np.matrix(_y).T
    
    return c
In [205]:
#eval svd error

actual,= plt.plot(np.linspace(0,1,N), 2+3*np.linspace(0,1,N))
plt.scatter(x, y)
axes = plt.gca()
axes.set_xlim([0,1])

coef = SVD(x, y)

print "svd:"
print "a = " + str(np.array(coef[0])[0][0])
print "b = " + str(np.array(coef[1])[0][0])
print ""

svd,= plt.plot(np.linspace(0,1,N), np.array(coef[0])[0][0]+np.array(coef[1])[0][0]*np.linspace(0,1,N))

#equation 12.34
a_error = stdev*pow(np.power(np.array(V[0,:])[0], 2).dot(np.power(s_inv, 2)),0.5)
b_error = stdev*pow(np.power(np.array(V[1,:])[0], 2).dot(np.power(s_inv, 2)),0.5)
print "12.34 error:"
print "error a = " + str(a_error) + "   error b = " + str(b_error)
print ""

import random
import math

#bootstrap

num_bootstrap = 100
boostrap_as = []
boostrap_bs = []
for i in range(num_bootstrap):

    bootstrap_x = [];
    bootstrap_y = [];
    for i in range(N):
        index = int(math.floor(random.random()*100))
        bootstrap_x.append(x[index])
        bootstrap_y.append(y[index])

    bootstrap_x = np.array(bootstrap_x)
    bootstrap_y = np.array(bootstrap_y)

    # actual,= plt.plot(np.linspace(0,1,N), 2+3*np.linspace(0,1,N))
    # plt.scatter(bootstrap_x, bootstrap_y)
    # axes = plt.gca()
    # axes.set_xlim([0,1])
    # plt.legend([actual], ['y = 2 + 3x'], loc=2)
    # plt.show()

    c = SVD(bootstrap_x, bootstrap_y)
    boostrap_as.append(np.array(c[0])[0][0])
    boostrap_bs.append(np.array(c[1])[0][0])

mean_a = np.array(boostrap_as).mean()
var_a = np.array(boostrap_as).var()
mean_b = np.array(boostrap_bs).mean()
var_b = np.array(boostrap_bs).var()

print "Bootstrap error:"
print "mean_a = " + str(mean_a)
print "mean_b = " + str(mean_b)
print "var a = " + str(var_a) + "   var b = " + str(var_b)
print ""

bootstrap,= plt.plot(np.linspace(0,1,N), mean_a+mean_b*np.linspace(0,1,N))


#many trials

num_trials = 100
aa = []
bb = []
for i in range(num_trials):

    _x = np.random.rand(N)
    _y = 2+3*_x
    for i in range(N):
        _y[i] += np.random.normal(0, stdev)
    
    # actual,= plt.plot(np.linspace(0,1,N), 2+3*np.linspace(0,1,N))
    # plt.scatter(bootstrap_x, bootstrap_y)
    # axes = plt.gca()
    # axes.set_xlim([0,1])
    # plt.legend([actual], ['y = 2 + 3x'], loc=2)
    # plt.show()

    c = SVD(_x, _y)
    aa.append(np.array(c[0])[0][0])
    bb.append(np.array(c[1])[0][0])

mean_a = np.array(aa).mean()
var_a = np.array(aa).var()
mean_b = np.array(bb).mean()
var_b = np.array(bb).var()

ensemble,= plt.plot(np.linspace(0,1,N), mean_a+mean_b*np.linspace(0,1,N))


print "Ensemble error:"
print "mean_a = " + str(mean_a)
print "mean_b = " + str(mean_b)
print "var a = " + str(var_a) + "   var b = " + str(var_b)


plt.legend([actual, svd, bootstrap, ensemble], ['y = 2 + 3x', 'Single SVD', "Bootstrap", "Ensemble"], loc=2)
plt.show()
svd:
a = 1.76274588475
b = 3.36748746285

12.34 error:
error a = 0.104751104519   error b = 0.176694975829

Bootstrap error:
mean_a = 1.7399152408
mean_b = 3.39458125174
var a = 0.00667950648253   var b = 0.0216714445469

Ensemble error:
mean_a = 1.98308245498
mean_b = 3.04201082781
var a = 0.0103967604105   var b = 0.0301651729061

Problem 12.2

In [207]:
N = 100
stdev = 0.1
x = np.random.rand(N)
y = np.sin(2+3*x)
for i in range(N):
    y[i] += np.random.normal(0, stdev)

actual,= plt.plot(np.linspace(0,1,N), np.sin(2+3*np.linspace(0,1,N)))
plt.scatter(x, y)
axes = plt.gca()
axes.set_xlim([0,1])
plt.legend([actual], ['y = sin(2 + 3x)'], loc=1)
plt.show()
In [214]:
actual,= plt.plot(np.linspace(0,1,N), np.sin(2+3*np.linspace(0,1,N)))


params = [1,1]

x_fit = np.linspace(0,1,N)

def evaluate():
    y_fit = np.sin(params[0]+params[1]*x_fit)
    fit,= plt.plot(x_fit, y_fit)

evaluate()

plt.scatter(x, y)
axes = plt.gca()
axes.set_xlim([0,1])
plt.legend([actual], ['y = sin(2 + 3x)'], loc=1)
plt.show()