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

<matplotlib.figure.Figure at 0x12d91e6d8>
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)

0.0369266215245
0.116188178639

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

<matplotlib.figure.Figure at 0x1200c0d30>
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()

<matplotlib.figure.Figure at 0x12b20e278>
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

In [777]:
def LevMarq(x,y,lam,a):
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)

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


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:
#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]:
0.90400592086568554
In [799]:
plt.plot(ahist[:,0])
plt.plot(ahist[:,1])

Out[799]:
[<matplotlib.lines.Line2D at 0x121027ba8>]
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]:
[<matplotlib.lines.Line2D at 0x12683df60>]
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):

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