problém odlehlých hodnot
M-odhady - pomocí vhodné váhové funkce je menší citlivost na vychýlené body než u nejm. čtverců ($d^2$)
minimalizujeme $\sum g(x_i,\theta)$, ve většině případů lze derivovat podle odhadu $\theta$: $\psi(x_i,\theta)=\partial g(x_i,\theta)/\partial \theta$ a váha se vyjádří jako
$$w(r_i)=\frac{\psi(r_i)}{r_i}$$ kde $r_i=x_i-\theta$
%matplotlib inline
from matplotlib import pyplot as pl
import matplotlib
matplotlib.rcParams['figure.figsize'] = [10, 5]
#mpl.use("pgf")
#mpl.use('GTK')
from numpy import r_,random
vytvoříme sadu dat s normálně rozdělenou nejistotou a příměsí rovnoměrně rozdělených vychýlených bodů
"správná" funkce je lineární se sklonem slp=-0.6
import numpy as np
frac=2.6
slp=-.6
x=r_[2:10:50j]
y=x*slp+random.normal(3,0.4,size=x.shape) #sigma=0.4
x2=np.concatenate([x,random.uniform(2,10,size=int(frac*x.size))])
y2=np.concatenate([y,random.uniform(y.min()-3,y.max()+3,size=int(frac*x.size))])
pl.plot(x2,y2,'.')
c=2 #tuning constant
gfunw=lambda z:c**2/2*(1-np.exp(-z**2/c**2)) #welsh
gfunh=lambda z:(abs(z)<c)*z**2/2+(abs(z)>=c)*(abs(c*z)-c**2/2) #huber
gfun=gfunh
#gfun=lambda z:abs(z)
wei=lambda p:gfun(y2-p[1]-p[0]*x2).sum()
from scipy import optimize as op
eslp=-0.5 # uvodni odhad sklonu
op.fmin(wei,[eslp,y2.mean()-eslp*x2.mean()])
#klasicke nejmensi ctverce
gfun=lambda z:z**2/2
wei=lambda p:gfun(y2-p[1]-p[0]*x2).sum()
op.fmin(wei,[eslp,y2.mean()-eslp*x2.mean()])
xn=r_[-5:5:0.1]
pl.plot(xn,gfunw(xn))
pl.plot(xn,gfunh(xn))
pl.legend(["Welsh","Huber"])
np.polyfit(x2,y2,1)
pl.ylabel("g-function")
pl.grid()
# jina volba tuning constant
c=1.
op.fmin(wei,[-0.5,y2.mean()-0.5*x2.mean()])
# nejmensi ctverce
gfunc=lambda z:z**2/2
gfun=gfunc
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
rep=[]
for j in range(len(y),len(y2),10):
x3,y3=x2[:j],y2[:j]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arep=np.array(rep)
arep
# median
gfund=lambda z:abs(z)
gfun=gfund
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
rep=[]
for j in range(len(y),len(y2),10):
x3,y3=x2[:j],y2[:j]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arepd=np.array(rep)
arepd
# Geman-McClure - vychazi podobne jako Welsch
gfunm=lambda z:z**2/2/(1+z**2)
gfun=gfunm
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
rep=[]
for j in range(len(y),len(y2),10):
x3,y3=x2[:j],y2[:j]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arepm=np.array(rep)
arepm
#Huberova funkce
gfun=gfunh
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for j in range(len(y),len(y2),10):
x3,y3=x2[:j],y2[:j]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
areph=np.array(rep)
areph
#Welshova funkce
gfun=gfunw
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for j in range(len(y),len(y2),10):
x3,y3=x2[:j],y2[:j]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arepw=np.array(rep)
arepw
#vysledky testu pro sklon
ix=np.arange(len(y),len(y2),10.)/len(y)
ix=1-1/ix
pl.plot(ix,arep[:,0],'s')
pl.plot(ix,arepw[:,0],'d')
pl.plot(ix,areph[:,0],'v')
pl.plot(ix,arepd[:,0],'o')
pl.legend(["quad","Welsh","Huber","absol."],loc=2)
pl.ylabel("slope")
pl.ylim(-0.7,0.1)
pl.grid()
pl.xlabel("frac. outlier.")
pl.xlim(-0.1,0.8)
pl.plot(ix,arep[:,1],'s')
pl.plot(ix,arepw[:,1],'d')
pl.plot(ix,areph[:,1],'v')
pl.plot(ix,arepd[:,1],'o')
pl.legend(["quad","Welsh","Huber"],loc=3)
pl.ylabel("absol.")
#pl.ylim(-0.7,-0.1)
pl.grid()
pl.xlabel("frac. outlier.")
pl.xlim(-0.1,0.8)
gfun=gfunh
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for c in np.r_[0.1:3:0.1]:
x3,y3=x2[:120],y2[:120]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arech=np.array(rep)
gfun=gfunw
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for c in np.r_[0.1:3:0.1]:
x3,y3=x2[:120],y2[:120]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
arecw=np.array(rep)
pl.plot(np.r_[0.1:3:0.1],arech[:,0],'gv')
pl.plot(np.r_[0.1:3:0.1],arecw[:,0],'yd')
pl.xlabel("tuning constant c")
pl.grid()
pl.legend(["Huber","Welsh"],loc=2)
co je bootstrap: vybereme bootfr bodů ze sady dat a nahradíme je nějakými (jinými) z původního vzorku
použitelné tehdy, pokud klasické metody (derivace minim. funkce podle parametru) selhávají
boofr=20
gsel=(len(x)+len(x2))//2
gfun=gfunw
rep=[]
for i in range(1000):
x3,y3=x2[:gsel],y2[:gsel].copy()
ifrom=np.random.choice(gsel,boofr)
ito=np.random.choice(gsel,boofr)
y3[ito]=y3[ifrom]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
brepw=np.array(rep)
brepw.std(0)
poměr dobrých a vychýlených dat 50:65
gfun=gfunc
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for i in range(1000):
x3,y3=x2[:gsel],y2[:gsel].copy()
ifrom=np.random.choice(gsel,boofr)
ito=np.random.choice(gsel,boofr)
y3[ito]=y3[ifrom]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
brepc=np.array(rep)
brepc.std(0)
gfun=gfunh
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for i in range(1000):
x3,y3=x2[:gsel],y2[:gsel].copy()
ifrom=np.random.choice(gsel,boofr)
ito=np.random.choice(gsel,boofr)
y3[ito]=y3[ifrom]
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
breph=np.array(rep)
breph.std(0)
# definujeme obecnou funkci
def bootme(gfunx, boofr=20,niter=1000,gsel=115):
rep=[]
wei=lambda p:gfunx(y3-p[1]-p[0]*x3).sum()
for i in range(niter):
x3,y3=x2[:gsel],y2[:gsel].copy()
ifrom=np.random.choice(gsel,boofr)
mask=np.ones(gsel,'bool')
mask[ifrom]=False
ito=np.random.choice(np.r_[:gsel][mask],boofr)
y3[ito]=y3[ifrom]
#gfun=gfunx
wei=lambda p:gfunx(y3-p[1]-p[0]*x3).sum()
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars),y3.std()])
return np.array(rep)
bootme(gfunh,20,500).std(0)
c=2
bootme(gfunw,20,500,115).std(0)
# Welsch
bootme(gfunw,20,500,115).std(0)
#zamenujeme 2x tolik bodu
bootme(gfunw,40,500,115).std(0)
zavisí na různých parametrech - bootstrap slouží jako relativní hodnota pro porovnání metod
# McClure vychazi nejlepe..
barr=bootme(gfunm,20,500)
barr.std(0),barr.mean(0)
#kolik ruznych prvku bude mit nahodny vyber s opakovanim
reprod=np.array([len(set(np.random.randint(gsel,size=gsel))) for i in range(200)])
reprod.mean(),reprod.min(),reprod.max(),reprod.std()
numerický experiment místo bootstrapu
gfun=gfunw
x1=x.copy()
c=1
rep=[]
for i in range(500):
y1=x1*slp+random.normal(3,0.4,size=x.shape)
x3=np.concatenate([x1,random.uniform(2,10,size=int(gsel-x.size))])
y3=np.concatenate([y1,random.uniform(y.min()-3,y.max()+3,size=int(gsel-x.size))])
pars=op.fmin(wei,[eslp,y3.mean()-eslp*x3.mean()],disp=0)
rep.append(list(pars)+[wei(pars)])
crepw=np.array(rep)
crepw.std(0)
rep=[]
for c in np.r_[0.6:3:0.2]:
barr=bootme(gfunw,20,500,115)
rep.append(np.r_[barr.std(0),barr.mean(0)])
bootdis=np.array(rep)
#linearni clen
bootdis[:,[4,0]]
#absolutni clen
bootdis[:,[5,1]]
gfun2=lambda x:abs(x)
wei2=lambda p:np.median(gfun2(y2-p[1]-p[0]*x2))
op.fmin(wei2,[-0.5,y.mean()-0.5*x.mean()])
Odstranění rozdělením vzorku (M rozdělíme na $m_1$ a $m_2$ o velikostech N)
$$E(M)=\theta + \frac{1}{2N}\beta + O(1/N^2)$$
$$E(m_{1/2})=\theta + \frac{1}{N}\beta + O(1/N^2)$$
$$2E(M)-\frac{E(m_1)+E(m_2)}{2} = \theta + O(1/N^2)$$
Zvýší se ovšem rozptyl členem řádu $1/N$
optimalizovaná varianta, rozptyl naroste jen člen řádu $1/N^2$ (porovnání s metodou bootstrapu)
$$\hat \alpha^* = N \hat \alpha - \frac{N-1}{N}\sum \alpha_j,$$ kde $\alpha_j=f(x_1,x_2,..x_{j-1},x_{j+1}...)$ je odhad s vynecháním j-tého vzorku.
ořez (trimming) - vynechání $n/2$ nejmenších a největších hodnot
winsorizace - nahrazení ořezaných hodnot minimem/maximem zbytku (kombinace průměru a midrange)
(location parameter)
posuzuje se asymptotická efektivita (poměr k MVB pro $N \to \infty$)
pro normálně rozdělená data s rostoucím $n$ efektivita klesá, pro rozdělení s větším vlivem okrajových hodnot (dvojexponenciála, Cauchy) roste
Optimum (minimax) pro $r=(N-n)/2N=0.23$
size=1000
def gendata(gaus=2,extr=10,exran=5,exshi=1):
data=random.normal(0,2,size)
if extr>0: data[::extr]+=random.uniform(-exran+exshi,exran+exshi,size//extr) #every extr-th point is an outlier
return data
data=gendata()
moms=[(data**i).mean() for i in [1,2]]
from numpy import sqrt
vari=moms[1]-moms[0]**2
[sum((data-moms[0])**2>i**2*vari) for i in [1,2,3]]
def bootest(nite=200,nsig=2,ntrim=0,extr=10,exran=10):
'''vraci prumer a smerod. odchylku z
- prumeru
- orezaneho prumeru na nsig nasobek smerod. odchylky
'''
rep=[]
for j in range(nite):
bootstr=gendata(extr=10,exran=exran)[random.randint(0,size,size)] #
moms=[bootstr.mean(),(bootstr**2).mean()] # vypocet momentu
vari=moms[1]-moms[0]**2 #odhad rozptylu (vychyleny?)
orez=bootstr[(bootstr-moms[0])**2<nsig**2*vari].mean()
if ntrim>0:
mtrim=np.sort(bootstr)[ntrim:-ntrim].mean()
rep.append([moms[0],orez,mtrim])
else:
rep.append([moms[0],orez])
arep=r_[rep]
return arep.mean(0),arep.std(0)
bootest(ntrim=100)
# pro 2 sigma a 10% orez
zoo=np.array([bootest(ntrim=100) for i in range(200)])
#dual exponent - trimming efficiency
q=np.r_[0.01:0.5:0.01]
D=2-(2-2*np.log(q)+np.log(q)**2)*q
pl.plot(q,1/D)
abs(zoo).mean(0)
pl.hist(abs(zoo[:,0,0]),r_[:0.03:20j],alpha=.5)
pl.hist(abs(zoo[:,0,1]),r_[:0.03:20j],alpha=.5)
# porovnani bez vychylenych bodu (extr=0)
zoo2=r_[[bootest(extr=0) for i in range(200)]]
abs(zoo2).mean(0)
zoo3=np.array([bootest(nsig=1,ntrim=200) for i in range(200)])
# pro 1 sigma a 20% orez
zoo3.mean(0)