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"