from numpy import*
import numpy as np

#Use an SVD to fit y = ax + b to the generated data set
N=100
M=2
seed = array([1,0.01])
x = zeros([N])
gamma = zeros([N])
for i in range (0,N-1):
	seed = ((8121*seed + 28411)%134456)
	x[i] = seed[0]/134456
	x[i+1] = seed[1]/134456
for i in range (0,N-1):
    seed = ((8121*seed + 28411)%134456)
    gamma[i] = 0.5*sqrt(-2*log(seed[0]/134456))*sin(seed[1])
    gamma[i+1] = 0.5*sqrt(-2*log(seed[0]/134456))*cos(seed[1])
y = 2+3*x + gamma

def fit(x,y):
	F = ones([N,M])
	F[:,M-1] = x
	u,w,v = np.linalg.svd(F, full_matrices = 0)
	W = diag(w)
	#values of a and b:
	A = dot(dot(dot(transpose(v),linalg.inv(W)),transpose(u)),y)
	return A, u, w, v
A,u,w,v = fit(x,y)
print "\na and b from SVD = \n" + str(A) + "\n"
	
#(a) Using equation (12.34) to evaluate the errors in a and b
errors = zeros([M])
for i in range(0,2):
	errors[i] = 0.25*sum(v[i,:]**2/w**2) #element-wise division
print "Errors in a and b using equation (12.34) = \n" + str(errors) + "\n"

#(b) Using bootstrapping (generating 100 data sets) to evaluate the errors in a and b
seed = 1
BS_x = zeros([N])
BS_y = zeros([N])
BS_A = zeros([100,M])
for i in range (0,100):
	for j in range (0,N):
		seed = ((8121*seed + 28411)%134456)	
		index = (int)(N*seed/(134456))
		BS_x[j] = x[index]
		BS_y[j] = y[index]
	A,u,w,v = fit(BS_x,BS_y)
	BS_A[i,:] = A;
M1 = array([sum(BS_A[:,0])/len(BS_A[:,0]),sum(BS_A[:,1])/len(BS_A[:,1])])
M2 = array([sum(BS_A[:,0]**2)/len(BS_A[:,0]),sum(BS_A[:,1]**2)/len(BS_A[:,1])])
BS_errors = M2-M1**2
print "Errors in a and b using bootstrapping = \n" + str(BS_errors) + "\n"
	
#(c) Using an ensemble of 100 independent data sets
ensemble_A = zeros([100,M])
seed = array([1,0.01])
for k in range(0,100):
	for i in range (0,N-1):
		seed = ((8121*seed + 28411)%134456)
		x[i] = seed[0]/134456
		x[i+1] = seed[1]/134456
	for i in range (0,N-1):
		seed = ((8121*seed + 28411)%134456)
		gamma[i] = 0.5*sqrt(-2*log(seed[0]/134456))*sin(seed[1])
		gamma[i+1] = 0.5*sqrt(-2*log(seed[0]/134456))*cos(seed[1])
	y = 2+3*x + gamma	
	A,u,w,v = fit(x,y)	
	ensemble_A[k,:] = A;
M1 = array([sum(ensemble_A[:,0])/len(ensemble_A[:,0]),sum(ensemble_A[:,1])/len(ensemble_A[:,1])])
M2 = array([sum(ensemble_A[:,0]**2)/len(ensemble_A[:,0]),sum(ensemble_A[:,1]**2)/len(ensemble_A[:,1])])
ensemble_errors = M2-M1**2
print "Errors in a and b using ensemble = \n" + str(ensemble_errors) + "\n"