#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()