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')
    gradient = np.zeros((2,1))
    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)
        grad = -8*(Y-sin)*cos
        gradient[0,0] = np.sum(grad)
        gradient[1,0] = np.sum(grad*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 gradient
        #print M
        #print np.linalg.pinv(M)
          
        dA = np.linalg.solve(-M,gradient)
        #print M, gradient, dA        
        errors.append(gradient)
        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))
        gradient = np.zeros((2,1))
        grad = -8*(Y-sin)*cos
        gradient[0,0] = np.sum(grad)
        gradient[1,0] = np.sum(grad*X)
        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
        return M,gradient,-np.linalg.pinv(M).dot(gradient)

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)[1][0] for a in np.linspace(0,4,100)],label='grad')
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 [ ]: