## SVD Line Fitting¶

### y = 2 + 3x + $\zeta$¶

In [2]:
import numpy as np
import matplotlib.pyplot as plt
import sympy

In [419]:
X = np.random.rand(100).reshape((100,1))

In [420]:
a_param=2
b_param=3

Y = a_param+b_param*X + np.random.normal(0,0.5,(100,1))

In [421]:
A = np.concatenate([np.ones((100,1)),X],axis=1)

A.shape

Out[421]:
(100, 2)
In [422]:
W=np.zeros((2,100))

U,s,Vt = np.linalg.svd(A)

W[0,0]=1/s[0];
W[1,1]=1/s[1];

In [423]:
a = Vt.T.dot(W).dot(U.T).dot(Y)[:2,:]
Y_est = A.dot(a)

a_est = a[0,0]
b_est = a[1,0]

In [425]:
plt.figure(figsize=(12, 8))
plt.plot(X,Y,'x')
plt.plot(X,Y_est, )
plt.show()


### Variance in a and b relative to variance in the data¶

In [426]:
a_bootstraps = np.zeros((2,100))
a_rands = np.zeros((2,100))

for i in range(100):
inds = np.floor(np.random.rand(100)*100).astype('int')
W=np.zeros((2,100))
U,s,Vt = np.linalg.svd(A[inds,:])
W[0,0]=1/s[0];
W[1,1]=1/s[1];
a_bootstraps[:2,i] = Vt.T.dot(W).dot(U.T).dot(Y[inds,:])[:2,0]

for i in range(100):
W=np.zeros((2,100))
X_rand = np.random.rand(100).reshape((100,1))
Y_rand = a_param+b_param*X_rand + np.random.normal(0,0.5,(100,1))
A_rand = np.concatenate([np.ones((100,1)),X_rand],axis=1)
U,s,Vt = np.linalg.svd(A_rand)
W[0,0]=1/s[0];
W[1,1]=1/s[1];
a_rands[:2,i] = Vt.T.dot(W).dot(U.T).dot(Y_rand)[:2,0]

print "SVD variance est:", np.sum(Vt.T**2/s**2, axis=1)/0.5**2
print "resample variance est:", np.var(a_bootstraps, axis=1)/0.5**2
print "new sample variance est:", np.var(a_rands, axis=1)/0.5**2

SVD variance est: [ 0.18555203  0.49379402]
resample variance est: [ 0.03657201  0.11683532]
new sample variance est: [ 0.03661118  0.1126607 ]


## Levenberg-Marquardt Fitting a Sinusoidal function¶

### y = sin(2 + 3x) + $\zeta$¶

In [6]:
import numpy as np
import matplotlib.pyplot as plt

X = np.random.rand(100)
Y = np.sin(2+3*X)+np.random.normal(0,0.1,100)

In [29]:
def Levenberg(X, Y, a=1.5, b=1.5, n=500, threshold=0.1, maxIters=1000, vary=False):
diff = float('inf')
M = np.zeros((2,2))
i = 0;
errors = [];

while (diff>threshold and i<maxIters):
oldDiff = diff
sin = np.sin(a+b*X)
cos = np.cos(a+b*X)

m = 4*(cos**2+(Y-sin)*sin)

M[0,0] = np.sum(m)*(1+n)
M[[0,1],[1,0]] = np.sum(m*X)
M[1,1] = np.sum(m*X**2)*(1+n)

diff = sum(((Y-sin)/0.5)**2)

if vary:
if diff>oldDiff:
n*=1.5
else:
n/=1.4

#print M
#print np.linalg.pinv(M)

a += dA[0]
b += dA[1]
i += 1
return a,b,errors


In [34]:
def test(X, Y, a, b):
sin = np.sin(a+b*X)
cos = np.cos(a+b*X)
M = np.zeros((2,2))
m = 4*(cos**2+(Y-sin)*sin)
diff = sum(((Y-sin)/0.5)**2)
M[0,0] = np.sum(m)*20
M[[0,1],[1,0]] = np.sum(m*X)
M[1,1] = np.sum(m*X**2)*20

plt.figure(figsize=(8, 6))
plt.title('M, gradient, and dA for different values of a')

plt.plot(np.linspace(0,4,100),[test(X,Y,a,1)[0][0,0] for a in np.linspace(0,4,100)],label='M')
plt.plot(np.linspace(0,4,100),[test(X,Y,a,1)[2][0] for a in np.linspace(0,4,100)],label='dA')
plt.plot(np.linspace(0,4,100),[0 for x in range(100)])
plt.ylim(-10,100)
plt.legend()
plt.show()


### Able to converge to a,b if starting close. Unable to converge from 1,1¶

In [37]:
from ipywidgets import *
%matplotlib notebook

def show(a_start, b_start):
a,b,errors = Levenberg(X,Y,a=a_start,b=b_start, n=50)

plt.figure(figsize=(8, 6))
plt.plot(X,Y,'x')
X2 = np.linspace(0,5,100)
print 'start:', a_start, b_start
print 'end:', a, b
plt.title('lambda=50')
plt.plot(X2,np.sin(a+b*X2),label='levenberg params')
plt.plot(X2,np.sin(a_start+b_start*X2),ls='--',label='start params')
plt.legend()
plt.show()

interact(show,a_start=(0,4,0.2), b_start=(0,4,0.2))

start: 1.2 1.8
end: [ 1.96233974] [ 3.08398603]

In [36]:
from ipywidgets import *
%matplotlib notebook

def show(a_start, b_start):
a,b,errors = Levenberg(X,Y,a=a_start,b=b_start, n=50, vary=True)
plt.figure(figsize=(8, 6))
plt.plot(X,Y,'x')
plt.title('lambda=50 at start but varies')
X2 = np.linspace(0,5,100)
print 'start:', a_start, b_start
print 'end:', a, b
plt.plot(X2,np.sin(a+b*X2),label='levenberg params')
plt.plot(X2,np.sin(a_start+b_start*X2),ls='--',label='start params')
plt.legend()
plt.show()

interact(show,a_start=(0,4,0.2), b_start=(0,4,0.2))

[ 1.96878379] [ 3.06957638]

In [ ]: