modely v zobecněném pojetí

  • kernelové vyhlazování
In [1]:
%matplotlib inline
from numpy import *
nen=30
x=r_[:10:1j*nen]+random.normal(size=nen)*0.1
y=exp(-(x-3)**2+0.1*(x-3)**3)+0.4*exp(-(x-6)**2)
y+=random.normal(size=nen)*0.1
In [3]:
y.max()
Out[3]:
0.95491800174285413
In [4]:
from matplotlib import pyplot as plt
plt.plot(x,y,'d')
fun=lambda x:exp(-(x-3)**2+0.1*(x-3)**3)+0.4*exp(-(x-6)**2)
plt.plot(r_[:10:.1],fun(r_[:10:.1]))
plt.grid()

klouzavý průměr odpovídá konstantnímu kernelu

průměrování obecně

  • snižuje výšku píků
  • zvyšuje šířku
  • zachovává momenty - plochu pod píkem a těžiště (nultý a první mom.)
In [5]:
k=5
xf=r_[:10:100j]
vf=[]
vf8=[]
for f in xf:
    ford=argsort(abs(f-x))
    vf.append(y[ford[:5]].mean())
    vf8.append(y[ford[:8]].mean())
plt.plot(xf,vf)
plt.plot(xf,vf8)
plt.legend(["5 sousedu","8 sousedu"])
plt.ylim(-0.4,1)
plt.grid()

nyní zkusíme "vážené" průměrování - kernel klesající se vzdáleností od středu

In [6]:
k=5
xf=r_[:10:100j]
kf=[]
kf1=[]
kf2=[]
for f in xf:
    dist=abs(f-x)
    ford=argsort(dist)
    norm=exp(-dist[ford[:10]]**2/0.2)
    norm/=norm.sum()
    norm1=exp(-dist[ford[:10]]**2/0.4)
    norm1/=norm1.sum()
    norm2=(1-dist[ford[:10]])
    norm2[norm2<0]=0
    norm2/=norm2.sum()
    kf.append((y[ford[:10]]*norm).sum())
    kf1.append((y[ford[:10]]*norm1).sum())
    kf2.append((y[ford[:10]]*norm2).sum())
    #vf8.append(y[ford[:8]].mean())
plt.plot(xf,kf)
plt.plot(xf,kf1)
plt.plot(xf,kf2)
plt.legend(["gauss 0.1","gauss 0.2","triangl"])
plt.ylim(-0.4,1)
plt.grid()

průměrování po skupinách

In [7]:
ymid=y.reshape(10,3).mean(1)
ystd=y.reshape(10,3).std(1)/sqrt(3)
xmid=x.reshape(10,3).mean(1)
xstd=x.reshape(10,3).std(1)/sqrt(3)

ymid2=y.reshape(6,5).mean(1)
ystd2=y.reshape(6,5).std(1)/sqrt(5)
xmid2=x.reshape(6,5).mean(1)
xstd2=x.reshape(6,5).std(1)/sqrt(5)

plt.errorbar(xmid,ymid,ystd,xstd)
plt.errorbar(xmid2,ymid2,ystd2,xstd2)
Out[7]:
<Container object of 3 artists>

fitovani analytickymi (parametrickymi) funkcemi

In [8]:
fun2=lambda p:p[0]+p[1]*exp(-(xc-3)**2)+p[2]*exp(-(xc-6)**2)
fun3=lambda p:p[0]+p[1]*exp(-(xc-3)**2+0.1*(xc-3)**3)+p[2]*exp(-(xc-6)**2)
fun3x=lambda p:p[0]+p[1]*exp(-(xc-3)**2+p[3]*(xc-3)**3)+p[2]*exp(-(xc-6)**2)
mat2=array([ones(x.shape),exp(-(x-3)**2),exp(-(x-6)**2)])
mat3=array([ones(x.shape),exp(-(x-3)**2+0.1*(x-3)**3),exp(-(x-6)**2)])
  • model 2: 2 symetrické píky
  • model 3: správná asymetrie 1. píku
In [9]:
dcov2=linalg.inv(mat2.dot(mat2.T)) #hessian
dcov3=linalg.inv(mat3.dot(mat3.T))
dcov2.dot(mat2.dot(y)),dcov3.dot(mat3.dot(y))
Out[9]:
(array([-0.0304616 ,  1.01125507,  0.35176036]),
 array([-0.03198407,  1.01478937,  0.3436464 ]))
In [10]:
pars2,pars3=dcov2.dot(mat2.dot(y)),dcov3.dot(mat3.dot(y))
plt.plot(x,y,'d')
xc=r_[:10:0.1]
plt.plot(xc,fun2(pars2))
plt.plot(xc,fun3(pars3))
plt.grid()

porovnání s nelineární optimalizací

In [11]:
from scipy import optimize
xc=x
minpars2=optimize.fmin(lambda p:((y-fun2(p))**2).sum(),[0,1,1])
minpars3=optimize.fmin(lambda p:((y-fun3(p))**2).sum(),[0,1,1])
pars2-minpars2,pars3-minpars3

Optimization terminated successfully.
         Current function value: 0.432061
         Iterations: 111
         Function evaluations: 188
Optimization terminated successfully.
         Current function value: 0.411611
         Iterations: 112
         Function evaluations: 201
Out[11]:
(array([  3.94483606e-06,   1.69843852e-05,  -2.65554972e-05]),
 array([ -1.42914194e-05,   3.77224705e-05,   3.87506706e-05]))
In [12]:
minpars3x=optimize.fmin(lambda p:((y-fun3x(p))**2).sum(),[0,1,1,.3])
minpars3x

Optimization terminated successfully.
         Current function value: 0.409029
         Iterations: 185
         Function evaluations: 312
Out[12]:
array([-0.03324669,  1.0156079 ,  0.34018385,  0.12671995])
In [13]:
#skutecne reziduum
((y-fun(x))**2).sum()
Out[13]:
0.46711036260139466
In [14]:
mat1=array([exp(-(x-3)**2),exp(-(x-6)**2)])
dcov1=linalg.inv(mat1.dot(mat1.T))
pars1=dcov1.dot(mat1.dot(y))
pars1,((y-mat1.T.dot(pars1))**2).sum()
Out[14]:
(array([ 0.96666419,  0.30749253]), 0.44646423997524776)
In [15]:
sqrt(dcov2.diagonal()),sqrt(dcov1.diagonal())
Out[15]:
(array([ 0.25382123,  0.66359584,  0.65529695]),
 array([ 0.54982497,  0.54162276]))