Fitovani polynomem 3. stupne

41 merenych bodu (interval (-8,4))

In [1]:
import numpy as np
a0,b0,c0,d0=3.,4.,2.,1.
nlev=0.2
x=np.r_[-8:4:41j][1:]
ydata=lambda x:a0+b0*x+c0*x**2+d0*x**3+np.random.normal(0,nlev,size=x.shape)
amat=np.array([np.ones_like(x),x,x**2,x**3])
hess=amat.dot(amat.T)
cov=np.linalg.inv(hess)
def ahat(y):
    pars=cov.dot(amat.dot(y))
    s0=((y-amat.T.dot(pars))**2).sum()
    return list(pars)+[s0]
pars=np.array([ahat(ydata(x)) for i in range(1000)]) #1000 iteraci
pars.mean(0)
Out[1]:
array([2.99865974, 3.99990248, 2.00031737, 1.00005234, 1.42725174])
In [2]:
#korelacni matice
errs=np.sqrt(cov.diagonal())
corx=cov/errs.reshape(4,1)/errs.reshape(1,4)
corx
Out[2]:
array([[ 1.        ,  0.24100027, -0.74565286, -0.59701386],
       [ 0.24100027,  1.        , -0.2450737 , -0.6111973 ],
       [-0.74565286, -0.2450737 ,  1.        ,  0.87744901],
       [-0.59701386, -0.6111973 ,  0.87744901,  1.        ]])
In [3]:
#porovnani odhadu nejistoty s rozptylem iterovanych experimentu
errs/pars.std(0)[:4]
Out[3]:
array([5.16287384, 4.9953337 , 4.96591069, 4.97018885])
In [4]:
# generujeme sit parametru 11x11x11x11
npt=11
zbas=slice(-0.05,0.05,1j*npt)
zand=np.mgrid[zbas,zbas,zbas,zbas]
y=ydata(x) # jeden pokus generovani dat
zpars=ahat(y)
alpars=np.array([zpars[i]*(1+zand[i]).ravel() for i in range(4)]).T
#alchis=np.array([(difarr(c)**2).sum() for c in alpars]).reshape(npt,npt,npt,npt)
alpars.min(0),alpars.max(0)
Out[4]:
(array([2.90355944, 3.81451048, 1.89852719, 0.94940221]),
 array([3.20919727, 4.2160379 , 2.09837215, 1.04933929]))
In [5]:
alchi=np.array([((p.dot(amat)-y)**2).sum() for p in alpars])-zpars[-1]
%matplotlib inline
from matplotlib import pyplot as pl
gchi=alchi.reshape(npt,npt,npt,npt)
pl.imshow(gchi.mean(0).mean(2))
pl.colorbar()
alchi.min()
Out[5]:
0.0
In [7]:
pl.imshow(gchi[6,:,:,6])
pl.colorbar()
Out[7]:
<matplotlib.colorbar.Colorbar at 0x7fdd99a7b160>
In [17]:
pl.plot(pars[:,1],pars[:,3],'.')
pl.xlabel("param. 2")
pl.ylabel("param. 4");
In [18]:
jl=np.r_[:4].astype(int)
allids=np.concatenate([np.r_[[np.ones(3-i)*i,jl[jl>i]]].T for i in jl]).astype(int)
zandz=zand+np.array(zpars[:4]).reshape(4,1,1,1,1)
modmat=list((zand**2).reshape(4,npt**4))
#modmat=list(alpars.T**2)
for ids in allids:
    modmat.append((zand[ids[0]]*zand[ids[1]]).ravel())
#    modmat.append((alpars[:,ids[0]]*alpars[:,ids[1]]).ravel())
modmat=np.array(modmat)
modmat.shape
Out[18]:
(10, 14641)
In [9]:
allids
Out[9]:
array([[0, 1],
       [0, 2],
       [0, 3],
       [1, 2],
       [1, 3],
       [2, 3]])
In [11]:
pl.plot(gchi.mean(0).mean(2).mean(1))
pl.plot(gchi.mean(0).mean(2).mean(0))
pl.legend(["param. 2","param. 3"])
Out[11]:
<matplotlib.legend.Legend at 0x7fdd99b7aa90>
In [9]:
%matplotlib inline
from matplotlib import pyplot as pl
gchi=alchi.reshape(npt,npt,npt,npt)
pl.imshow(gchi.mean(0).mean(2))
pl.colorbar()
Out[9]:
<matplotlib.colorbar.Colorbar instance at 0x7f0baa1cdb00>
In [11]:
hessx=modmat.dot(modmat.T)
covx=np.linalg.inv(hessx)
parx=covx.dot(modmat.dot(alchi.ravel()))
parx
Out[11]:
array([  3.61828373e+02,   9.82218264e+03,   8.22552400e+04,
         8.81619298e+05,  -1.77658206e+03,   7.39967629e+03,
        -1.75293338e+04,  -4.64320727e+04,   1.64879320e+05,
        -5.06274268e+05])
In [12]:
pl.plot(modmat[0].reshape(npt,npt*npt*npt).mean(1))
parx[:4]*errs**2
Out[12]:
array([ 27.05565478,  79.41414751,  77.82600622,  20.84971242])
In [12]:
modmat.shape
Out[12]:
(10, 14641)
In [13]:
dermat=np.ones((4,4))
dermat[jl,jl]=parx[:4]
for i in range(len(allids)):
    dermat[allids[i][0],allids[i][1]]=parx[i+4]/2.
    dermat[allids[i][1],allids[i][0]]=parx[i+4]/2.
#(covmat*100).astype(int)
covmat=np.linalg.inv(dermat)
covmat
Out[13]:
array([[  8.26633119e-03,   4.93647131e-04,  -1.04523942e-03,
         -2.64097411e-04],
       [  4.93647131e-04,   5.07557662e-04,  -8.51259737e-05,
         -6.69957811e-05],
       [ -1.04523942e-03,  -8.51259737e-05,   2.37708871e-04,
          6.58215212e-05],
       [ -2.64097411e-04,  -6.69957811e-05,   6.58215212e-05,
          2.36726384e-05]])
In [14]:
cov/covmat
Out[14]:
array([[  9.04570933,  12.00393283,   6.00038622,   3.00612803],
       [ 12.00393283,  15.92958586,   7.96269596,   3.9892238 ],
       [  6.00038622,   7.96269596,   3.98029977,   1.99408676],
       [  3.00612803,   3.9892238 ,   1.99408676,   0.99901571]])
In [15]:
gvals=np.r_[a0,b0,c0,d0]
cov/covmat/gvals.reshape(1,4)/gvals.reshape(4,1)
Out[15]:
array([[ 1.00507881,  1.00032774,  1.00006437,  1.00204268],
       [ 1.00032774,  0.99559912,  0.995337  ,  0.99730595],
       [ 1.00006437,  0.995337  ,  0.99507494,  0.99704338],
       [ 1.00204268,  0.99730595,  0.99704338,  0.99901571]])
In [16]:
dif = lambda p:((p.dot(amat)-y)**2).sum()
qpars=np.r_[zpars[:4]]
In [62]:
def estmat(msize=1000,npar=4,eran=0.05,loud=0):
    pts=np.random.uniform(-eran,eran,size=npar*msize).reshape(npar,msize)
    mpts=(1+pts)*qpars[:,np.newaxis]
    chi0=dif(qpars)
    allchi=[dif(b)-chi0 for b in mpts.T]
    if loud==-2: return mpts
    if loud==-1: return allchi
    if loud: print(min(allchi),max(allchi))
    jl=np.r_[:npar].astype(int)
    #allids=np.concatenate([np.r_[[np.ones(npar-1-i)*i,jl[jl>i]]].T for i in jl]).astype(int)
    #zandz=zand+np.array(zpars[:npar]).reshape(npar,1,1,1,1)
    zand=pts
    modmat=list((zand**2).reshape(npar,msize))
    for ids in allids:
        modmat.append((zand[ids[0]]*zand[ids[1]]).ravel())
    modmat=np.array(modmat)
    if loud==2: return modmat
    hessx=modmat.dot(modmat.T)
    covx=np.linalg.inv(hessx)
    parx=covx.dot(modmat.dot(allchi))
    if loud==3: return parx
    if loud==4: return parx,((allchi-parx.dot(modmat))**2).sum()
    dermat=np.ones((npar,npar))
    dermat[jl,jl]=parx[:npar]
    for i in range(len(allids)):
        dermat[allids[i][0],allids[i][1]]=parx[i+npar]/2.
        dermat[allids[i][1],allids[i][0]]=parx[i+npar]/2.
    #(covmat*100).astype(int)
    covmat=np.linalg.inv(dermat)
    return covmat
covmat/estmat(1500)#*gvals.reshape(1,4)*gvals.reshape(4,1)/cov
Out[62]:
array([[ 103.51627902, -271.64058916,  764.56321251, -607.65233211],
       [-271.64058916,    5.84937218,   17.86686092,    7.74893644],
       [ 764.56321251,   17.86686092,    3.4909028 ,    4.6053145 ],
       [-607.65233211,    7.74893644,    4.6053145 ,    1.73623605]])
In [64]:
zpars,ahat(y)
Out[64]:
([3.005352490361588,
  3.9742849500067337,
  2.0050883182106958,
  1.0011450875657566,
  1.426621020196964],
 [2.1198148916137143,
  3.7670860204047187,
  2.054784362583149,
  1.0100828479018329,
  112.6141083493456])
In [63]:
dif = lambda p:((p.dot(amat)-y)**2).sum()
qpars=np.r_[zpars[:4]]
covmat/estmat(100)
Out[63]:
array([[-7.96942255, -1.87956948, -5.42592591, -3.77110656],
       [-1.87956948,  1.84371841, -6.00843666,  2.7681983 ],
       [-5.42592591, -6.00843666,  3.88883096,  3.61490132],
       [-3.77110656,  2.7681983 ,  3.61490132,  2.53579661]])
In [22]:
covmat/estmat(50)
Out[22]:
array([[ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.],
       [ 1.,  1.,  1.,  1.]])
In [81]:
rep=[]
for j in range(50):
    y=ydata(x) # jeden pokus generovani dat
    zpars=ahat(y)
    dif = lambda p:((p.dot(amat)-y)**2).sum()
    qpars=np.r_[zpars[:4]]
    newmat=estmat(100)
    nerrs=np.sqrt(newmat.diagonal())
    newcor=newmat/nerrs.reshape(1,4)/nerrs.reshape(4,1)
    iarr=np.r_[[[0,1],[0,2],[0,3],[1,2],[1,3],[2,3]]]
    ia,ib=iarr.T
    rep.append(np.r_[nerrs*qpars,newcor[ia,ib]])
arep=np.array(rep).T
arep.mean(1)
Out[81]:
array([ 0.27344987,  0.08991765,  0.03075959,  0.00486306,  0.24100027,
       -0.74565286, -0.59701386, -0.2450737 , -0.6111973 ,  0.87744901])
In [43]:
arep.std(1)
Out[43]:
array([  4.48762562e-12,   1.23489979e-13,   2.80863157e-13,
         2.85309323e-14,   5.02857392e-12,   5.47499482e-12,
         6.45572896e-12,   4.25989458e-12,   1.56371798e-12,
         1.56053946e-12])
In [44]:
arep.mean(1)
Out[44]:
array([ 0.27344987,  0.08991765,  0.03075959,  0.00486306,  0.24100027,
       -0.74565286, -0.59701386, -0.2450737 , -0.6111973 ,  0.87744901])

now some robust methods

In [46]:
a=1.51
rho=lambda t:t**2*(abs(t)<a)+a*(abs(t)>a)*abs(t)
phi=lambda t:t*(abs(t)<a)+a*(abs(t)>a)*np.sign(t)
dphi=lambda t: (abs(t)<a)
In [ ]:
dif = lambda p:((p.dot(amat)-y)**2).sum()
In [54]:
msel=np.random.uniform(size=y.shape[0])<0.1
y[msel]+=np.random.uniform(-10,10,sum(msel))
In [65]:
from scipy import optimize as op
dif=lambda p:rho(p.dot(amat)-y).sum()
op.fmin(dif,qpars)

Optimization terminated successfully.
         Current function value: 25.243248
         Iterations: 87
         Function evaluations: 148
Out[65]:
array([ 2.93039505,  3.95720143,  2.00940828,  1.00190771])
In [82]:
dif2=lambda p:rho(p.dot(amat)-y).sum()
npt=11
zoo=np.ones((npt,4))*qpars
zoo[:,2]*=(1+np.r_[-0.015:0.015:1j*npt])
pl.plot((1+np.r_[-0.015:0.015:1j*npt]),[dif2(z) for z in zoo])
Out[82]:
[<matplotlib.lines.Line2D at 0x7f0ba6382750>]
In [95]:
dif1=lambda p:((p.dot(amat)-y)**2).sum()
y=ydata(x)
qpars0=np.r_[ahat(y)[:4]]
#y[msel]+=np.random.uniform(-5,5,sum(msel))
dif=dif2
qpars=qpars0
estmat(50,eran=0.02)/covmat
Out[95]:
array([[ 0.00542251, -0.00050412, -0.00176551,  0.00137742],
       [-0.00050412,  0.13324191,  0.04300988,  0.11240633],
       [-0.00176551,  0.04300988,  0.12470111,  0.10133734],
       [ 0.00137742,  0.11240633,  0.10133734,  0.26886502]])
In [107]:
dif=dif2
y[msel]+=np.random.uniform(-5,5,sum(msel))
allchi=np.array(estmat(50,eran=0.02,loud=-1))
allchi.min(),allchi.max()
Out[107]:
(0.81002111469451776, 113.74600424651453)