x = np.linspace(0,2,50)
y = x + x**2 + np.random.randn(len(x))

V = matrix(np.vstack([ones(x.shape),x,x**2]).T)
Q = np.matrix(np.eye(V.shape[0]))
i,j =np.diag_indices_from(Q)
Q[i[:20],j[:20]]=100
R = np.matrix(np.diag([1,1,13]))*0

Pv = V*inv(V.T*Q*V+R)*V.T*Q

p=np.polyfit(x,y,2)

fig, ax = subplots()
fig.set_size_inches(5,5)

ax.plot(x,y,'o',label='data',alpha=0.3)
ax.plot(x,np.dot(Pv,y).flat,label='projection')
ax.plot(x,np.polyval(p,x),'-',label='polyfit',alpha=0.3)
ax.grid()
ax.legend(loc=0)

x = np.random.rand(50)*2
x.sort() # sort for plotting
       
y = x + x**2 + np.random.randn(len(x))

V = matrix(np.vstack([ones(x.shape),x,x**2]).T)
Pv = V*inv(V.T*V+eye(3)*4)*V.T

p=np.polyfit(x,y,2)

fig, ax = subplots()
fig.set_size_inches(5,5)

ax.plot(x,y,'o',label='data',alpha=0.3)
ax.plot(x,np.dot(Pv,y).flat,'-s',label='projection')
#ax.plot(x,np.polyval(p,x),'s',label='polyfit')
ax.grid()
ax.legend(loc=0);

import numpy as np
import matplotlib.pyplot as plt
from matplotlib.widgets import Slider 

fig, ax = plt.subplots()
fig.set_size_inches(5,5)

plt.subplots_adjust(left=0.25, bottom=0.30,right=1)

x = np.linspace(0,2,50)
y = x + x**2 + np.sin(2*np.pi*x) + np.random.randn(len(x))

V = np.matrix(np.vstack([np.ones(x.shape),x,x**2]).T)
Q = np.matrix(np.eye(V.shape[0]))
i,j =np.diag_indices_from(Q)
#Q[i[:20],j[:20]]=100
R = np.matrix(np.diag([1,1,1]))*0

Pv = V*np.linalg.inv(V.T*Q*V+R)*V.T*Q
p=np.polyfit(x,y,2)

ax.plot(x,y,'o',label='data',alpha=0.3)
ax.set_title('err^2=%3.2f'%(np.linalg.norm(y)**2 - np.linalg.norm(np.dot(Pv,y))**2 ))
ls_line,=ax.plot(x,np.dot(Pv,y).flat,label='projection')
ax.plot(x,np.polyval(p,x),'-',label='polyfit',alpha=0.3)
ax.legend(loc=0)
ax.grid()

axw0 = plt.axes([0.25, 0.2, 0.65, 0.03] )
sw0  = Slider(axw0, 'w0', 0.1, 30.0, valinit = 1)
axw1 = plt.axes([0.25, 0.15, 0.65, 0.03] )
sw1  = Slider(axw1, 'w1', 0.1, 30.0, valinit = 1)
axw2 = plt.axes([0.25, 0.10, 0.65, 0.03] )
sw2  = Slider(axw2, 'w2', 0.1, 30.0, valinit = 1)

def update_w0(val):
   w = sw0.val
   R[0,0]=w
   Pv = V*np.linalg.inv(V.T*Q*V+R)*V.T*Q
   ls_line.set_ydata(np.dot(Pv,y).flat)
   ax.set_title('err=%3.2f'%(np.linalg.norm(y)**2 - np.linalg.norm(np.dot(Pv,y))**2 ))
   plt.draw()
sw0.on_changed(update_w0)

def update_w1(val):
   w = sw1.val
   R[1,1]=w
   Pv = V*np.linalg.inv(V.T*Q*V+R)*V.T*Q
   ls_line.set_ydata(np.dot(Pv,y).flat)
   ax.set_title('err=%3.2f'%(np.linalg.norm(y)**2 - np.linalg.norm(np.dot(Pv,y))**2 ))
   plt.draw()
sw1.on_changed(update_w1)

def update_w2(val):
   w = sw2.val
   R[2,2]=w
   Pv = V*np.linalg.inv(V.T*Q*V+R)*V.T*Q
   ls_line.set_ydata(np.dot(Pv,y).flat)
   ax.set_title('err=%3.2f'%(np.linalg.norm(y)**2 - np.linalg.norm(np.dot(Pv,y))**2 ))
   plt.draw()
sw2.on_changed(update_w2)

plt.show()