fixed 9_2

week9/week9_11_2.py

@@ -19,20 +19,28 @@ def generate_data(size, sigma):
 
 # Levenberg-Marquardt routine to fit y = sin(a+bx) to this data set starting from a = b = 1
 def levenberg_step(x, y, a, b, l):
-    # Calculate the Jacobian
-    J = np.array([np.ones(len(x)), x]).T
-    # Calculate the Hessian
-    H = np.dot(J.T, J)
-    # Calculate the gradient
-    g = np.dot(J.T, y - np.sin(a + b*x))
-    # update = (H + l*I)^-1 * g (where I is the identity matrix np.eye(2)
-    coefficients = np.dot(np.linalg.inv(H + l*np.eye(2)), g)
-    return coefficients
+    # every noisy y - every y from our fit:
+    error = (y - np.sin(a + b*x))**2
+    error_sum = np.sum(error)
+    J = np.zeros(2)
+    J[0] = -2 * np.sum((y - np.sin(a + b*x)) * np.cos(a + b * x))    # from wolfram
+    J[1] = -2 * np.sum((y - np.sin(a + b*x)) * np.cos(a + b * x) * x)
+
+    M = np.zeros((2, 2))
+    M[0, 0] = 0.5 * (1 + l) * 2 * np.sum(np.cos(a + b * x) ** 2)    # from wolfram
+    M[1, 1] = 0.5 * (1 + l) * 2 * np.sum((x * np.cos(a + b * x)) ** 2)
+
+    M[0, 1] = 0.5 * 2 * np.sum(x * (np.cos(a + b * x) ** 2)) 
+    M[1, 0] = M[0, 1] # it's a symmetric matrix
+
+    M_inv = np.linalg.inv(M)
+    step = -M_inv.dot(J)
+    return step
 
 a = 1
 b = 1
-l = 0.1
-steps = 100
+l = 1
+steps = 50
 
 size = 100
 sigma = 0.1
@@ -45,11 +53,12 @@ for i in range(steps):
     coefficients = levenberg_step(x, y, a, b, l)
     A.append(a)
     B.append(b)
-    a += coefficients[0]
+    a += coefficients[0]  # we are taking little relative moves of a and b, not returning absolute positions
     b += coefficients[1]
     # print("a: %.4f, b: %.4f" % (a, b))
 
-'''
+print(a,b)
+# '''
 # Here we plot the data and the fit line
 fig, (ax1, ax2) = plt.subplots(2, 1)
 # slim plot
@@ -64,8 +73,9 @@ ax1.set_title('Levenberg-Marquardt Fit\nsteps = %.2f, sigma = %.2f' % (steps, si
 ax2.scatter(A, B, c=range(steps))
 ax2.set_title('Converging a and b values')
 plt.show()
-'''
+# '''
 
+'''
 # Here we animate the fit converging over the scatter plot
 my_dpi = 96
 fig = plt.figure(figsize=(800/my_dpi, 400/my_dpi), dpi=my_dpi)
@@ -83,9 +93,10 @@ def animate(w):
     y = B[w]
     fit[0].set_data(x_sorted, np.sin(x + y*x_sorted))
 
-anim = animation.FuncAnimation(fig, animate, frames=int(20))
+anim = animation.FuncAnimation(fig, animate, frames=int(50))
 # plt.show()
 
 path = r"c://Users/david/Desktop/NMM/week9/img/animation.gif" 
 # save as a gif
 anim.save(path, writer='imagemagick', fps=8, dpi=my_dpi)
+'''
\ No newline at end of file