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$
from matplotlib import pyplot as pl
#mpl.use("pgf")
#mpl.use('GTK')
%matplotlib inline
from numpy import r_,random
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)
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
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
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)
c=1.
op.fmin(wei,[-0.5,y2.mean()-0.5*x2.mean()])
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
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
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
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.legend(["quad","Welsh","Huber"],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.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)
vybereme bootfr bodů ze sady dat a nahradíme je nějakými (jinými) z původního vzorku
boofr=20
j=(len(x)+len(x2))/2
gfun=gfunw
rep=[]
for i in range(1000):
x3,y3=x2[:j],y2[:j].copy()
ifrom=np.random.choice(j,boofr)
ito=np.random.choice(j,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)
gfun=gfunc
rep=[]
wei=lambda p:gfun(y3-p[1]-p[0]*x3).sum()
for i in range(1000):
x3,y3=x2[:j],y2[:j].copy()
ifrom=np.random.choice(j,boofr)
ito=np.random.choice(j,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[:j],y2[:j].copy()
ifrom=np.random.choice(j,boofr)
ito=np.random.choice(j,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)
gfun=gfunh
def bootme(boofr=20,niter=1000,gsel=115):
rep=[]
wei=lambda p:gfun(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]
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(20,500).std(0)
gfun=gfunw
bootme(20,500,115).std(0)
bootme(40,500,115).std(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()
brepw.mean(0), breph.mean(0)
gfun=gfunw
x1=x.copy()
rep=[]
for i in range(100):
y1=x1*slp+random.normal(3,0.4,size=x.shape)
x3=np.concatenate([x1,random.uniform(2,10,size=int(j-x.size))])
y3=np.concatenate([y1,random.uniform(y.min()-3,y.max()+3,size=int(j-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)
j
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):
data=random.normal(0,2,size)
if extr>0: data[::extr]+=random.uniform(-exran,exran,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,extr=10,exran=10):
rep=[]
for j in range(200):
bootstr=gendata(extr=10,exran=exran)[random.randint(0,size,size)]
moms=[(bootstr**i).mean() for i in [1,2]]
vari=moms[1]-moms[0]**2
rep.append([moms[0],bootstr[(bootstr-moms[0])**2<nsig**2*vari].mean()])
arep=r_[rep]
return arep.mean(0),arep.std(0)
zoo=[bootest() for i in range(200)]
#dual exponent - trimming efficiency
import numpy as np
%matplotlib inline
from matplotlib import pyplot as pl
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[:,0,:]).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)
zoo2=r_[[bootest(extr=0) for i in range(200)]]
abs(zoo2).mean(0)