In [753]:
%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import collections, colors, transforms
import numpy as np
import seaborn as sns
from mpl_toolkits.mplot3d import Axes3D
from matplotlib import cm
from matplotlib.ticker import LinearLocator, FormatStrFormatter
import numpy as np
import sympy
from matplotlib import animation, rc
from IPython.display import HTML
import matplotlib.mlab as mlab
from scipy import stats
In [ ]:
#Generate 100 points x uniformly distributed between 0 and 1, and let y = 2+3x+ζ,
#where ζ is a Gaussian random variable with a standard deviation of 0.5. Use an
#SVD to fit y = a + bx to this data set, finding a and b. Evaluate the errors in a
#and b
#(a) With equation (12.34)
#(b) By bootstrap sampling to generate 100 data sets
#(c) From fitting an ensemble of 100 independent data sets
In [735]:
### GENERATE DATA
x = np.random.rand(100)
g = np.random.normal(0,0.5,100)
y =2+3*x+g
### FIT MODEL
def SVDfit(x,y):
A = np.ones((100,2))
A[:,1] = x
U, s, V = np.linalg.svd(A,full_matrices = False)
S = np.diag(s)
temp = (U.T).dot(y)
temp = (np.linalg.inv(S)).dot(temp)
C = (V.T).dot(temp)
return(C)
C = SVDfit(x,y)
T = np.ones((2,2))
T[:,1] = np.arange(0.,2.)
test = T.dot(C)
fig = plt.figure()
fig=plt.figure(figsize=(4, 4), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
ax.scatter(x,y, s = 10,color = 'k')
ax.plot(np.arange(0.,2.),test,'r-',linewidth = 2,label='model')
ax.plot(np.arange(0.,2.),(3*np.arange(0.,2.))+2,'b--',linewidth = 1,label='ground truth')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
plt.show()
In [638]:
### GET FIT ERROR USING EQUATION
error1 = np.sum((V[0,:]**2)/((s)**2))
error2 = np.sum((V[1,:]**2)/((s)**2))
print(error1)
print(error2)
In [736]:
### GET FIT ERROR USING BOOTSTRAPPING
### BOOTSTRAP DATA
def bootstrap(x,y):
X = np.ones((100,100))
Y = np.ones((100,100))
C = np.ones((100,2))
for i in range(0,100):
idx = np.random.randint(0,100,100)
X[:,i] = x[idx]
Y[:,i] = y[idx]
C[i,:] = SVDfit(X[:,i],Y[:,i])
return(X,Y,C)
X,Y,C = bootstrap(x,y)
fig = plt.figure()
fig=plt.figure(figsize=(4, 4), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
for i in range(0,99):
T = np.ones((2,2))
T[:,1] = np.arange(0.,2.)
test = T.dot(C[i,:])
ax.plot(np.arange(0.,2.),test,'r-',linewidth = 0.5)
T = np.ones((2,2))
T[:,1] = np.arange(0.,2.)
test = T.dot(C[99,:])
ax.plot(np.arange(0.,2.),test,'r-',linewidth = 0.5,label='bootstrapped models')
ax.scatter(x,y, s = 10,color = 'k')
ax.plot(np.arange(0.,2.),(3*np.arange(0.,2.))+2,'b--',linewidth = 1,label='ground truth')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
plt.show()
C1 = C
In [734]:
fig = plt.figure()
fig=plt.figure(figsize=(4, 4), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
C = np.ones((100,2))
for i in range(0,100):
x = np.random.rand(100)
g = np.random.normal(0,0.5,100)
y =2+3*x+g
C[i,:] = SVDfit(x,y)
for i in range(0,99):
T = np.ones((2,2))
T[:,1] = np.arange(0.,2.)
test = T.dot(C[i,:])
ax.plot(np.arange(0.,2.),test,'r-',linewidth = 0.5)
C2 = C
T = np.ones((2,2))
T[:,1] = np.arange(0.,2.)
test = T.dot(C[99,:])
ax.plot(np.arange(0.,2.),test,'r-',linewidth = 0.5,label='models')
ax.plot(np.arange(0.,2.),(3*np.arange(0.,2.))+2,'b--',linewidth = 1,label='ground truth')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
plt.show()
In [750]:
fig = plt.figure()
#fig=plt.figure(figsize=(4, 4), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
ax.scatter(np.zeros(100),C1[:,0],s=5,color = 'r')
ax.scatter(np.ones(100),C1[:,1],s=5,color = 'r',label='bootstrapped fits')
ax.scatter(2*np.ones(100),C2[:,0],s=5,color = 'b')
ax.scatter(3*np.ones(100),C2[:,1],s=5,color = 'b',label='independent fits')
ax.plot([-1,4],[2,2],'--k')
ax.plot([-1,4],[3,3],'--k',label='ground truth')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
plt.show()
In [756]:
fig = plt.figure()
#fig=plt.figure(figsize=(4, 4), dpi= 180, facecolor='w', edgecolor='k')
ax = fig.gca()
ax.plot(0,error1,'ok')
ax.plot(1,error2,'ok') #,label='ground truth')
ax.plot(2,np.var(C1[:,0]),'or')
ax.plot(3,np.var(C1[:,1]),'or')
ax.plot(4,np.var(C2[:,0]),'ob')
ax.plot(5,np.var(C2[:,1]),'ob')
handles, labels = ax.get_legend_handles_labels()
ax.legend(handles, labels)
plt.show()
In [250]:
#Generate 100 points x uniformly distributed between 0 and 1, and let y = sin(2 +
#3x) + ζ, where ζ is a Gaussian random variable with a standard deviation of 0.1.
#Write a Levenberg-Marquardt routine to fit y = sin(a+bx) to this data set starting
#from a = b = 1, and investigate the convergence for both fixed and adaptively
#adjusted λ values.
In [777]:
def LevMarq(x,y,lam,a):
grad_chi = [0,0]
M = np.zeros((2,2))
dya0 = np.cos(a[0] + a[1]*x)
dya1 = x*np.cos(a[0] + a[1]*x)
d2ya0a1 = -x*np.sin(a[0] + a[1]*x)
d2ya02 = -np.sin(a[0] + a[1]*x)
d2ya12 = -(x**2)*np.sin(a[0] + a[1]*x)
res = y - np.sin(a[0]+a[1]*x)
grad_chi[0]= -2* np.sum((res/np.var(y))*dya0)
grad_chi[1]= -2* np.sum((res/np.var(y))*dya1)
M[0,0] = (1+lam) * np.sum((1/np.var(y))*(dya0*dya0))# - res*d2ya02))
M[1,1] = (1+lam) * np.sum((1/np.var(y))*(dya1*dya1))# - res*d2ya12))
M[0,1] = M[1,0] = np.sum((1/np.var(y))*(dya0*dya1))# - res*d2ya0a1))
#da = -(np.linalg.inv(M)).dot(grad_chi)
return(M,grad_chi)
In [797]:
x = np.random.rand(100)
g = np.random.normal(0,0.1,100)
y = np.sin(2 + 3*x) + g
flag = True
eta = 1
lam = 10
a = [1,1]
error = [np.sum((y - np.sin(a[0]+a[1]*x))**2)]
ahist = a
dahist = [0,0]
lamhist = lam
j = 0
while flag and j < 10000:
M,grad_chi = LevMarq(x,y,lam,a)
a += -dot(np.linalg.inv(M),grad_chi)
#a = a-da
dahist = np.vstack((dahist,da))
ahist = np.vstack((ahist, a))
lamhist = np.vstack((lamhist, lam))
error = np.append(error,[np.sum((y - np.sin(a[0]+a[1]*x))**2)],axis = 0)
#if error is decreasing decrease lambda
if (error[-1]) - (error[-2]) < 0:
lam+= 1
else:
lam-= 1
j+=1
if error[-1] < 0.5:
flag = False
In [798]:
plt.plot(error)
error[-1]
Out[798]:
In [799]:
plt.plot(ahist[:,0])
plt.plot(ahist[:,1])
Out[799]:
In [801]:
plt.plot(x,y,'k.')
plt.plot(np.linspace(0,1,1000),np.sin(2+3*np.linspace(0,1,1000)),'b--')
plt.plot(np.linspace(0,1,1000),np.sin(a[0]+a[1]*np.linspace(0,1,1000)),'r-')
Out[801]:
In [802]:
fig = plt.figure()
ax = fig.add_subplot(111, autoscale_on=False, xlim=(0, 1), ylim=(-1, 1))
ax.grid()
line, = ax.plot([], [], 'r-', lw=2)
line1, = ax.plot([], [], 'b-', lw=2)
#scatter, = ax.scatter([],[], 'k.')
time_template = 'time = %.1fs'
time_text = ax.text(0.05, 0.9, '', transform=ax.transAxes)
def init():
line.set_data([], [])
line1.set_data([],[])
#scatter.set_data([],[])
return line, line1#, scatter
def animate(i):
thisx = [np.linspace(0,1,1000)]
thisy = [np.sin(ahist[i,0]+ahist[i,1]*np.linspace(0,1,1000))]
#scatter.set_data(x,y)
line.set_data(thisx, thisy)
line1.set_data(np.linspace(0,1,1000),np.sin(2+3*np.linspace(0,1,1000)))
return line, line1#, scatter
anim = animation.FuncAnimation(fig, animate, 100,interval=50, blit=True, init_func=init)
plt.close()
HTML(anim.to_html5_video())
Out[802]:
In [473]:
def LevMarq(x,y,lam,a):
grad_chi[0] = -2*np.sum(((y - np.sin(a[0]+a[1]*x))/np.var(y))*(np.cos(a[0]+a[1]*x)))
grad_chi[1] = -2*np.sum(((y - np.sin(a[0]+a[1]*x))/np.var(y))*(x*np.cos(a[0]+a[1]*x)))
M[0,0] = (1+lam)* np.sum((1/np.var(y))*(
((np.cos(a[0]+a[1]*x))**2)-
(y - np.sin(a[0]+a[1]*x))*
(-np.sin(a[0]+a[1]*x))))
M[1,1] = (1+lam)* np.sum((1/np.var(y))*(
((x*(np.cos(a[0]+a[1]*x)))**2)-
(y - np.sin(a[0]+a[1]*x))*
(-x**2*np.sin(a[0]+a[1]*x))))
M[1,0] = np.sum((1/np.var(y))*(((np.cos(a[0]+a[1]*x)))*(x*np.cos(a[0]+a[1]*x))-(y - np.sin(a[0]+a[1]*x))*
(-x*np.sin(a[0]+a[1]*x))))
M[0,1] = np.sum((1/np.var(y))*(((np.cos(a[0]+a[1]*x)))*(x*np.cos(a[0]+a[1]*x))-(y - np.sin(a[0]+a[1]*x))*
(-x*np.sin(a[0]+a[1]*x))))
da = -(np.linalg.inv(M)).dot(grad_chi)
return(da,lam)