Metoda najmniejszych kwadratów

Aproksymacja wielomianowa

Zastosowanie algebry liniowej do prostego przypadku interpolacji funkcji wielomianami $x^n$.
Sage: Wykorzystanie technik operacji na macierzach i danych.

Problem:

Mamy dane punkty $(x_1,y_1),\dots,(x_n,y_n)$. Szukamy wielomianu stopnia $n-1$ przechodzącego przez te punkty. $c_i$ to współczynniki wielomianu



In [116]:

    
var('x')
xs=[]
ys=[]
cs=[]
n_pkts = 4
for i in range(1,n_pkts+1):
    xs.append(var('x%d' % i))
    ys.append(var('y%d' % i))
    cs.append(var('c%d' % i))
X=vector([x^i for i in range(n_pkts)]    )
A=(matrix([ X.subs(x==xx) for xx in xs] ))
b=(vector(ys).column())
c=(vector(cs).column())
# table([["$A$","$\cdot$","$c$","$= b$"],[A,"$\cdot$",c,b]],header_row=True)



In [117]:

    
show(A)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrrr} 1 & x_{1} & x_{1}^{2} & x_{1}^{3} \\ 1 & x_{2} & x_{2}^{2} & x_{2}^{3} \\ 1 & x_{3} & x_{3}^{2} & x_{3}^{3} \\ 1 & x_{4} & x_{4}^{2} & x_{4}^{3} \end{array}\right)$

Bierzemy konkretne punkty:



In [118]:

    
points=[[1,1],[2,2],[3,5]]
points









    Out[118]:





[[1, 1], [2, 2], [3, 5]]



In [119]:

    
table(  [["x","y"]]+points, header_row=True)
pkt_plt=point(points,size=30)
show(pkt_plt,figsize=4)



In [ ]:



In [120]:

    
npoints=len(points)
print npoints



In [121]:

    
[[x_^i for i in range(npoints)] for x_,y_ in points]









    Out[121]:





[[1, 1, 1], [1, 2, 4], [1, 3, 9]]



In [22]:

    
A = matrix( [[x_^i for i in range(npoints)] for x_,y_ in points])
show(A)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 1 & 1 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \end{array}\right)$



In [23]:

    
b = vector( [y_ for x_,y_ in points] )
b









    Out[23]:





(1, 2, 5)



In [ ]:

Rozwiązujemy ze względu na współczynniki $c_i$:



In [24]:

    
c = A.solve_right(b)
c









    Out[24]:





(2, -2, 1)

sprawdźmy:



In [25]:

    
A*c==b









    Out[25]:





True



In [ ]:



In [ ]:

    
[["$A$","$\cdot$","$c$","$= b$"]



In [32]:

    
show(A,c.column(),"=",b.column())









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 1 & 1 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \end{array}\right) \left(\begin{array}{r} 2 \\ -2 \\ 1 \end{array}\right) \verb|=| \left(\begin{array}{r} 1 \\ 2 \\ 5 \end{array}\right)$



In [33]:

    
## Możemy utworzyć  wielomian mnożąc skalarnie wektor z elementami bazy funckji ze współczynnikami



In [34]:

    
var('x')
X = vector([x^i for i in range(npoints)])



In [35]:

    
wielomian= c.dot_product(X)
show(wielomian)
plot(wielomian,(-.5,33.5),color='red',figsize=4)+pkt_plt
print(wielomian)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}x^{2} - 2 \, x + 2$ 






    



x^2 - 2*x + 2

Weźmy inny zestaw danych



In [36]:

    
n=12
l=[ ( 2*pi.n()/(n-1)*i   ,   sin(2*pi.n()/(n-1)*i)  ) for i in range(n)]


A = matrix( [[x_^i for i in range(n)] for x_,y_ in l])
b = vector( [y_ for x_,y_ in l] )
var('x')
X = vector([x^i for i in range(n)])
c = A.solve_right(b)
wielomian= c.dot_product(X)
pp = plot(wielomian,(x,-4,12),ymin=-2,ymax=2)
pp += plot(sin(x),(x,-4,4*pi),color='red')
pp += point(l,color='green',size=70)
pp.show()

W SageMath mamy instniejącą implementacje zagadnienia Lagrange'a.



In [115]:

    
var('x')
n=12
l=[ ( 2*pi.n()/(n-1)*i   ,   sin(2*pi.n()/(n-1)*i)  ) for i in range(n)]
R = PolynomialRing(RR, 'x')
L = R.lagrange_polynomial(l)
pp=L.plot(-4,12,ymin=-2,ymax=2)+plot(sin(x),(x,0,2*pi),color='red')+point(l,color='green',size=70)
pp.show()

Należy uważać na fakt, że x jest generatorem pierścienia wielomianów nie zmienną symboliczą:



In [118]:

    
var('y')
R.<x> = PolynomialRing(RR)
L = R.lagrange_polynomial(l)
print type(L(y))
print type(L)
print type(L(x))
print type(x)









    



<type 'sage.symbolic.expression.Expression'>
<type 'sage.rings.polynomial.polynomial_real_mpfr_dense.PolynomialRealDense'>
<type 'sage.rings.polynomial.polynomial_real_mpfr_dense.PolynomialRealDense'>
<type 'sage.rings.polynomial.polynomial_real_mpfr_dense.PolynomialRealDense'>



In [121]:

    
var('x')
f(x)=exp(-x)*sin(x)
f(1.1)









    Out[121]:





0.296657159743355



In [ ]:

interact

link

var('xx')
f(xx)=sin(xx)
@interact
def sinpoly(n=slider(range(2,20),default=4),a=slider(srange(0.2,5,0.2),default=1) ) :

    l=[ (a*2*pi.n()/(n-1)*i,exp(-0.2*a*2*pi.n()/(n-1)*i)*sin(a*2*pi.n()/(n-1)*i)) for i in range(n)]
    R = PolynomialRing(RR, 'x')
    L = R.lagrange_polynomial(l)
    plt=L.plot(-4,a*2*pi+4,ymin=-1.2,ymax=1.2)+plot(exp(-0.2*x)*sin(x),(x,0,2*pi*a),color='red')+point(l,color='green',figsize=(6,3),size=70)
    plt.show()

Inne bazy

Bazą nie muszą być wielomiany. Konieczne jest jednak by był to zestaw funkcji liniowo niezależnych.



In [181]:

    
points = [[1,1],[2,2],[3,5]]



In [182]:

    
x_lst = [x_ for x_,y_ in points]
y_lst = [y_ for x_,y_ in points]



In [183]:

    
X = vector([x^i for i in range(npoints)])

Spróbuj na przykład:

X = vector([sin(x*i) for i in range(npoints)])
    X[0] = 1

albo:

X = vector([exp(x*i) for i in range(npoints)])



In [184]:

    
X









    Out[184]:





(1, x, x^2)

Zapiszmy tworzenie macierzy $A$ w sposób niezależny od tego jaką formuła będzie w wektorze X:



In [185]:

    
A = matrix( [X.subs({x:x_}) for x_ in x_lst] ) 
b = vector(y_lst)



In [186]:

    
A









    Out[186]:





[1 1 1]
[1 2 4]
[1 3 9]



In [187]:

    
c = A.solve_right(b)



In [188]:

    
plot(X.dot_product(c),(x,-2,3)) + point(points,figsize=4,color='red')









    Out[188]:



In [ ]:

Metoda najmniejszych kwadratów

Cel:

Zrozumienie podstaw metody najmniejszych kwadratów
Sage: umiejętność operacji na macierzach i danych.

Zamiast:

$$ A x = b$$

rozwiązujemy

$$ A^T A x = A^T b.$$

Inaczej mówiąc, niech błąd:

$$ r=b-A x$$

leży w lewym jądrze A:

$$A^T r =A^T( b-Ax) = A^T b-A^TAx = 0.$$

$$ A^T A x = A^T b.$$

$$ \frac{\partial}{\partial x_k} (A_{ij} x_j - b_i) (A_{il} x_l - b_i) = 0$$

$$ A_{ij} \delta_{jk} (A_{il} x_l - b_i) + A_{il} \delta_{lk} (A_{ij} x_j - b_i) = $$ $$ A^TAx - A^T b + A^TAx - A^T b $$



In [179]:

    
var('x y t')
points=[[0,1],[2,0],[3,4],[5,7]]



In [180]:

    
N=len(points)

X=vector([1,x,x^2])
M=len(X)

var('x')
xs=[]
ys=[]
cs=[]
for i in range(1,N+1):
    xs.append(var('x%d' % i))
    ys.append(var('y%d' % i))
for i in range(1,M+1):    
    cs.append(var('c%d' % i))
    

A=(matrix([ X.subs(x==xx) for xx in xs] ))
b=(vector(ys).column())
c=(vector(cs).column())


ATA=A.transpose()*A
ATb=A.transpose()*b

#table([["$A$","$\cdot$","$c$","$= b$"],
       
show(A,c,"=",b)
show(ATA,c,"=",ATb)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & x_{1} & x_{1}^{2} \\ 1 & x_{2} & x_{2}^{2} \\ 1 & x_{3} & x_{3}^{2} \\ 1 & x_{4} & x_{4}^{2} \end{array}\right) \left(\begin{array}{r} c_{1} \\ c_{2} \\ c_{3} \end{array}\right) \verb|=| \left(\begin{array}{r} y_{1} \\ y_{2} \\ y_{3} \\ y_{4} \end{array}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 4 & x_{1} + x_{2} + x_{3} + x_{4} & x_{1}^{2} + x_{2}^{2} + x_{3}^{2} + x_{4}^{2} \\ x_{1} + x_{2} + x_{3} + x_{4} & x_{1}^{2} + x_{2}^{2} + x_{3}^{2} + x_{4}^{2} & x_{1}^{3} + x_{2}^{3} + x_{3}^{3} + x_{4}^{3} \\ x_{1}^{2} + x_{2}^{2} + x_{3}^{2} + x_{4}^{2} & x_{1}^{3} + x_{2}^{3} + x_{3}^{3} + x_{4}^{3} & x_{1}^{4} + x_{2}^{4} + x_{3}^{4} + x_{4}^{4} \end{array}\right) \left(\begin{array}{r} c_{1} \\ c_{2} \\ c_{3} \end{array}\right) \verb|=| \left(\begin{array}{r} y_{1} + y_{2} + y_{3} + y_{4} \\ x_{1} y_{1} + x_{2} y_{2} + x_{3} y_{3} + x_{4} y_{4} \\ x_{1}^{2} y_{1} + x_{2}^{2} y_{2} + x_{3}^{2} y_{3} + x_{4}^{2} y_{4} \end{array}\right)$



In [187]:

    
pkt_plt=point(points,size=30)
show(pkt_plt,figsize=4)



In [188]:

    
A = matrix( [X.subs(x==x_[0]) for x_ in points] ) 
b = vector( [x_[1] for x_ in points] ) 
show(A,b.column())









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 0 & 0 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \\ 1 & 5 & 25 \end{array}\right) \left(\begin{array}{r} 1 \\ 0 \\ 4 \\ 7 \end{array}\right)$



In [ ]:



In [189]:

    
ATA=A.transpose()*A
ATb=A.transpose()*b



In [191]:

    
show(A,c,"=",b.column())
show(ATA,c,"=",ATb.column())









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 1 & 0 & 0 \\ 1 & 2 & 4 \\ 1 & 3 & 9 \\ 1 & 5 & 25 \end{array}\right) \left(\begin{array}{r} c_{1} \\ c_{2} \\ c_{3} \end{array}\right) \verb|=| \left(\begin{array}{r} 1 \\ 0 \\ 4 \\ 7 \end{array}\right)$ 






    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left(\begin{array}{rrr} 4 & 10 & 38 \\ 10 & 38 & 160 \\ 38 & 160 & 722 \end{array}\right) \left(\begin{array}{r} c_{1} \\ c_{2} \\ c_{3} \end{array}\right) \verb|=| \left(\begin{array}{r} 12 \\ 47 \\ 211 \end{array}\right)$



In [192]:

    
c=ATA\ATb



In [193]:

    
c









    Out[193]:





(19/26, -14/39, 1/3)



In [194]:

    
wielomian= c.dot_product( X )
show(wielomian)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{3} \, x^{2} - \frac{14}{39} \, x + \frac{19}{26}$



In [195]:

    
pkt_plt+plot(wielomian, (x,0,5),color='red',thickness=2,figsize=5)









    Out[195]:

Praktyczne stosowanie interpolacji i aproksymacji w Sage

Zagadnienie Lagrange'a w Sage:



In [200]:

    
pts=[[1,1],[2,2],[3,4]]
R = PolynomialRing(QQ, 'x')    ### x - bedzie generatorem wielomianów nad ciałem liczb wymiernych, 
                               ### R - będzie objektem reprezentującym to ciało
L = R.lagrange_polynomial(pts) ### interpolacja Lagrange'a jest zaimplementowana w R 
show(L(x).expand() )









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, x^{2} - \frac{1}{2} \, x + 1$

Przykład: zjawisko Rungego

Rozważmy funkcję:



In [202]:

    
f(x)= 1/(1+25*x^2)
plot(f,(-1,1), figsize=4)









    Out[202]:

Jej pochodne przyjmują duże wartości w $0$:



In [203]:

    
[ diff(1/(1+25*x^2),x,i).subs({x:0}) for i in range(1,7) ]









    Out[203]:





[0, -50, 0, 15000, 0, -11250000]

Zjawisko Rungego - interact

link

@interact

    def inter1(n_pkt=slider(range(1,20),default=7) ):
        x0  = 1
        f(x)= 1/(1+25*x^2)
        x_bounds=(-1,1)
        p   = plot(f,x_bounds, thickness=2)

        #n_pkt=10
        xvals1=[cos(i/(n_pkt-1)*pi.n()) for i in range(n_pkt)]
        xvals2=[-1+ (2*i/(n_pkt-1)).n() for i in range(n_pkt)]
        #xvals3=[ 2*random()-1 for i in range(n_pkt) ]
        plt=[]
        for xvals in [xvals2,xvals1]:

            pts=[(i,QQ(f(i))) for i in xvals]
            R = PolynomialRing(QQ, 'x')
            L = R.lagrange_polynomial(pts)
            L_plot = plot(L(x),x_bounds, color='green', thickness=2,figsize=4)
            plt.append(L_plot+p+point(pts,pointsize=40,color='red'))

        print "           Punkty jednorodnie położone","                           Punkty Chebysheva"
        map(show,plt)



In [ ]:



In [204]:

    
pts=[(1,1),(2,2),(3,4)]



In [205]:

    
R = PolynomialRing(QQ, 'x')
L = R.lagrange_polynomial(pts)
show(L)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} x^{2} - \frac{1}{2} x + 1$



In [206]:

    
DD=R.divided_difference(pts)
DD









    Out[206]:





[1, 1, 1/2]



In [207]:

    
xlst=map(lambda y:y[0],pts)
xlst









    Out[207]:





[1, 2, 3]



In [211]:

    
f=[1]
for xi in xlst[:-1]:
    f.append( f[-1]*(x-xi) )
show(f)









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left[1, x - 1, {\left(x - 1\right)} {\left(x - 2\right)}\right]$



In [213]:

    
show(zip(DD,f))









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\left[\left(1, 1\right), \left(1, x - 1\right), \left(\frac{1}{2}, {\left(x - 1\right)} {\left(x - 2\right)}\right)\right]$



In [214]:

    
sum([dd*f_ for dd,f_ in zip(DD,f)]).expand().show()









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, x^{2} - \frac{1}{2} \, x + 1$



In [215]:

    
DD=vector(DD)
f=vector(f)
f.dot_product(DD).expand().show()









    




 $\newcommand{\Bold}[1]{\mathbf{#1}}\frac{1}{2} \, x^{2} - \frac{1}{2} \, x + 1$



In [92]:









    Out[92]:



In [87]:









    Out[87]:



In [54]:









    Out[54]:

Praktyczne stosowanie regresji liniowej i nie tylko

find_fit - używa pakietu scipy.optimize



In [217]:

    
points









    Out[217]:





[[0, 1], [2, 0], [3, 4], [5, 7]]



In [218]:

    
var('a,b,c,x')
model(x) = c*x^1+b
fit = find_fit(points,model,solution_dict=True)
fit









    Out[218]:





{b: -0.26923076923353717, c: 1.3076923076929787}



In [ ]:



In [219]:

    
model.subs(fit)









    Out[219]:





x |--> 1.3076923076929787*x - 0.26923076923353717



In [220]:

    
point(points)+plot(model.subs(fit),(x,0,10),color='red',figsize=5)









    Out[220]:



In [64]:









    Out[64]:

Splajny (GSL)



In [214]:

    
s = spline(points)
show(point(points) + plot(s,1,3, hue=.8,figsize=5))
s.list()









    












    Out[214]:





[[0, 1], [2, 0], [3, 4], [5, 6]]



In [215]:

    
GSLpoints









    Out[215]:





[[0, 1], [2, 0], [3, 4], [5, 6]]



In [216]:

    
v = [(i + sin(i)/2, i+cos(i^2)) for i in range(10)]
s = spline(v)
show(point(v) + plot(s,0,9, hue=.8,figsize=5))
s.list()









    












    Out[216]:





[(0, 1),
 (1/2*sin(1) + 1, cos(1) + 1),
 (1/2*sin(2) + 2, cos(4) + 2),
 (1/2*sin(3) + 3, cos(9) + 3),
 (1/2*sin(4) + 4, cos(16) + 4),
 (1/2*sin(5) + 5, cos(25) + 5),
 (1/2*sin(6) + 6, cos(36) + 6),
 (1/2*sin(7) + 7, cos(49) + 7),
 (1/2*sin(8) + 8, cos(64) + 8),
 (1/2*sin(9) + 9, cos(81) + 9)]

Bezier



In [217]:

    
pts=[(0,0),(.15,.1),(0.7,1.6),(1,0)]
path = [pts]
curve = bezier_path(path, linestyle='dashed', rgbcolor='green',figsize=5)
curve+point(pts)









    Out[217]:



In [66]:









    Out[66]: