From 0a0cbc2fe239ff83f615f3952237f46033695689 Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Sun, 1 Jan 2012 23:09:40 +0000 Subject: [PATCH] Small updates --- pywafo/src/wafo/kdetools.py | 213 +++++++++++++++++++++++++++--------- pywafo/src/wafo/wafodata.py | 16 +-- 2 files changed, 170 insertions(+), 59 deletions(-) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 8180e4e..1be6a4a 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -417,7 +417,8 @@ class _KDE(object): else: titlestr = 'Kernel density estimate (%s)' % self.kernel.name kwds2 = dict(title=titlestr) - kwds2['plot_kwds'] = dict(plotflag=1) + + kwds2['plot_kwds'] = kwds.pop('plot_kwds',dict(plotflag=1)) kwds2.update(**kwds) args = self.args if self.d == 1: @@ -1030,7 +1031,7 @@ class KRegression(_KDE): """ def __init__(self, data, y, p=0, hs=None, kernel=None, alpha=0.0, xmin=None, xmax=None, inc=128, L2=None): - self.tkde = TKDE(data, hs, kernel, alpha, xmin, xmax, inc, L2) + self.tkde = TKDE(data, hs=hs, kernel=kernel, alpha=alpha, xmin=xmin,xmax=xmax, inc=inc, L2=L2) self.y = y self.p = p @@ -3492,25 +3493,16 @@ def _get_data(n=100, symmetric=False, loc1=1.1, scale1=0.6, scale2=1.0): x = xi[i] y = yi[i] return x, y, fun1 -def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj', plotlog=False): - x,y, fun1 = _get_data(n, symmetric) - kreg_demo3(x,y,fun1, hs=None, fun='hisj', plotlog=False) - -def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): - import scipy.stats as st - - alpha=0.1 - z0 = -_invnorm(alpha/2) - + +def _get_regression_smooting(x,y,fun='hste'): hs1 = Kernel('gauss', fun=fun).get_smoothing(x) - n = x.size - hx = np.median(np.abs(x-np.median(x)))/0.6745*(4.0/(3*n))**0.2 + #hx = np.median(np.abs(x-np.median(x)))/0.6745*(4.0/(3*n))**0.2 if (y==True).any(): hs2 = Kernel('gauss', fun=fun).get_smoothing(x[y==True]) - hy = np.median(np.abs(y-np.mean(y)))/0.6745*(4.0/(3*n))**0.2 + #hy = np.median(np.abs(y-np.mean(y)))/0.6745*(4.0/(3*n))**0.2 else: hs2 = 4*hs1 - hy = 4*hx + #hy = 4*hx #hy2 = Kernel('gauss', fun=fun).get_smoothing(y) #kernel = Kernel('gauss',fun=fun) @@ -3518,7 +3510,20 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): #hopt = (hs1+4*hs2)/5 #kernel.get_smoothing(x) #hopt = hs2 hopt = sqrt(hs1*hs2) - #hopt=sqrt(hx*hy); + return hopt, hs1, hs2 +def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj', plotlog=False): + x,y, fun1 = _get_data(n, symmetric) + kreg_demo3(x,y,fun1, hs=None, fun='hisj', plotlog=False) + +def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): + import scipy.stats as st + + alpha=0.1 + z0 = -_invnorm(alpha/2) + + + n = x.size + hopt, hs1, hs2 =_get_regression_smooting(x,y,fun='hste') if hs is None: hs = hopt @@ -3528,12 +3533,12 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): #reverse = np.exp xmin, xmax = x.min(), x.max() - ni = max(2*int((xmax-xmin)/hopt)+1,5) + ni = max(2*int((xmax-xmin)/hopt)+3,5) print(ni) print(xmin, xmax) - - xi = np.linspace(xmin-hopt,xmax+hopt, ni) - xiii = np.linspace(xmin-hopt,xmax+hopt, 4*ni+1) + sml = hopt + xi = np.linspace(xmin-sml,xmax+sml, ni) + xiii = np.linspace(xmin-sml,xmax+sml, 4*ni+1) c = gridcount(x, xi) if (y==True).any(): @@ -3547,7 +3552,7 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): print("fact=%g" % (fact)) kreg = KRegression(x, y, hs=hs*fact, p=0) fiii = kreg(xiii) - yiii = stineman_interp(xiii, x, y) + yiii = stineman_interp(xiii, xi, yi) fit = fun1(xiii).clip(max=1.0) df = np.diff(fiii) eerr = np.abs(yiii-fiii).std()+ 0.5*(df[:-1]*df[1:]<0).sum()/n @@ -3608,12 +3613,14 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): # Jeffreys intervall a=b=0.5 #st.beta.isf(alpha/2, x+a, n-x+b) - ab = 0.055 + ab = 0.07 #0.055 pi1 = pi #fun1(xiii) pup2 = np.where(pi==1, 1, st.beta.isf(alpha/2, ciii*pi1+ab, ciii*(1-pi1)+ab)) plo2 = np.where(pi==0, 0, st.beta.isf(1-alpha/2, ciii*pi1+ab, ciii*(1-pi1)+ab)) + averr = np.trapz(pup2-plo2, xiii)/(xiii[-1]-xiii[0]) + 0.5*(df[:-1]*df[1:]<0).sum()/n + #f2 = kreg_demo4(x, y, hs, hopt) # Wilson score den = 1+(z0**2./ciii); xc=(pi1+(z0**2)/(2*ciii))/den; @@ -3623,8 +3630,8 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): #pup = (pi + z0*np.sqrt(pi*(1-pi)/ciii)).clip(min=0,max=1) # dont use #plo = (pi - z0*np.sqrt(pi*(1-pi)/ciii)).clip(min=0,max=1) - fg.plot(label='KReg grid') - f.plot(label='KRegression') + fg.plot(label='KReg grid' ) + f.plot(label='KReg averr=%2.3f ' %(averr)) labtxt = '%d CI' % (int(100*(1-alpha))) plt.fill_between(xiii, pup, plo, alpha=0.20,color='r', linestyle='--', label=labtxt) plt.fill_between(xiii, pup2, plo2,alpha = 0.20, color='b', linestyle=':',label='%d CI2' % (int(100*(1-alpha)))) @@ -3638,6 +3645,131 @@ def kreg_demo3(x,y, fun1, hs=None, fun='hisj', plotlog=False): plt.setp(h,yscale='log') #plt.show() return hs1, hs2 + + +def check_kreg_demo3(): + import wafo.fig as fig + + plt.ion() + k = 0 + for i, n in enumerate([50, 100,300,600, 4000]): + x,y, fun1 = _get_data(n, symmetric=True,loc1=1.0, scale1=0.6, scale2=1.25) + k0 = k + + for j, fun in enumerate(['hste']): + hsmax, hs1, hs2 =_get_regression_smooting(x,y,fun=fun) + for hi in np.linspace(hsmax*0.25,hsmax,9): + plt.figure(k) + k +=1 + unused = kreg_demo3(x,y,fun1, hs=hi, fun=fun, plotlog=False) + + #kreg_demo2(n=n,symmetric=True,fun='hste', plotlog=False) + fig.tile(range(k0,k)) + plt.ioff() + plt.show() + +def check_kreg_demo4(): + import wafo.fig as fig + + plt.ion() + #test_docstrings() + #kde_demo2() + #kreg_demo1(fast=True) + #kde_gauss_demo() + #kreg_demo2(n=120,symmetric=True,fun='hste', plotlog=True) + k = 0 + for i, n in enumerate([100,300,600]): + x,y, fun1 = _get_data(n, symmetric=True,loc1=0.6, scale1=0.3, scale2=1.25) + k0 = k + hopt1, h1,h2 = _get_regression_smooting(x,y,fun='hns') + hopt2, h1,h2 = _get_regression_smooting(x,y,fun='hste') + hopt = sqrt(hopt1*hopt2) + for j, fun in enumerate(['hste', 'hisj', 'hns', 'hstt']): + hsmax, hs1, hs2 =_get_regression_smooting(x,y,fun=fun) + + fmax = kreg_demo4(x, y, hsmax+0.1, hopt) + for hi in np.linspace(hsmax*0.25,hsmax,25): + f = kreg_demo4(x, y, hi, hopt) + if f.prediction_error_avg<=fmax.prediction_error_avg: + fmax = f + plt.figure(k) + k +=1 + fmax.plot() + xi = fmax.args[::4] + c = gridcount(x, xi) + if (y==True).any(): + c0 = gridcount(x[y==True],xi) + else: + c0 = np.zeros(xi.shape) + yi = np.where(c==0, 0, c0/c) + plt.plot(xi, yi, 'b.', xi, fun1(xi),'r') + + #kreg_demo2(n=n,symmetric=True,fun='hste', plotlog=False) + fig.tile(range(k0,k)) + plt.ioff() + plt.show() + +def kreg_demo4(x,y, hs, hopt): + import scipy.stats as st + + n = x.size + alpha=0.1 + z0 = -_invnorm(alpha/2) + + xmin, xmax = x.min(), x.max() + ni = max(2*int((xmax-xmin)/hopt)+3,5) + + sml = hopt + xi = np.linspace(xmin-sml,xmax+sml, ni) + xiii = np.linspace(xmin-sml,xmax+sml, 4*ni+1) + + + from wafo.interpolate import stineman_interp + + kreg = KRegression(x, y, hs=hs, p=0) + f = kreg(xiii, output='plotobj', plot_kwds = dict(plotflag=7)) + + + dx = xi[1]-xi[0] + ciiii = kreg.tkde.eval_grid_fast(xiii)*dx* x.size + ckreg = KDE(x,hs=hs) + ciii = ckreg.eval_grid_fast(xiii)*dx* x.size #n*(1+symmetric) + + + pi = f.data + + # Jeffreys intervall a=b=0.5 + #st.beta.isf(alpha/2, x+a, n-x+b) + ab = 0.07 #0.055 + pi1 = pi #fun1(xiii) + pup = np.where(pi1==1, 1, st.beta.isf(alpha/2, ciii*pi1+ab, ciii*(1-pi1)+ab)) + plo = np.where(pi1==0, 0, st.beta.isf(1-alpha/2, ciii*pi1+ab, ciii*(1-pi1)+ab)) + + + # Wilson score + den = 1+(z0**2./ciii); + xc=(pi1+(z0**2)/(2*ciii))/den; + halfwidth=(z0*sqrt((pi1*(1-pi1)/ciii)+(z0**2/(4*(ciii**2)))))/den + plo2 = (xc-halfwidth).clip(min=0) # wilson score + pup2 = (xc+halfwidth).clip(max=1.0) # wilson score + + f.dataCI = np.vstack((plo,pup)).T + f.prediction_error_avg = np.trapz(pup-plo, xiii)/(xiii[-1]-xiii[0]) + fiii = f.data + + c = gridcount(x, xi) + if (y==True).any(): + c0 = gridcount(x[y==True],xi) + else: + c0 = np.zeros(xi.shape) + yi = np.where(c==0, 0, c0/c) + yiii = stineman_interp(xiii, xi, yi) + df = np.diff(fiii) + eerr = np.abs(yiii-fiii).std()+ 0.5*(df[:-1]*df[1:]<0).sum()/n + f.labels.title='perr=%1.3f,eerr=%1.3f, n=%d, hs=%1.3f' % (f.prediction_error_avg,eerr,n,hs) + + return f + def kde_gauss_demo(n=50): ''' KDEDEMO Demonstrate the KDEgauss @@ -3697,35 +3829,14 @@ def test_docstrings(): doctest.testmod() if __name__ == '__main__': - import wafo.fig as fig + #check_kreg_demo3() + check_kreg_demo4() + + #test_smoothn_2d() + #test_smoothn_cardioid() - plt.ion() #test_docstrings() #kde_demo2() #kreg_demo1(fast=True) #kde_gauss_demo() #kreg_demo2(n=120,symmetric=True,fun='hste', plotlog=True) - k = 0 - for i, n in enumerate([50, 100,300,600, 4000]): - x,y, fun1 = _get_data(n, symmetric=True,loc1=0.5, scale1=0.3, scale2=.75) - k0 = k - - for j, fun in enumerate(['hste']): - hs1 = Kernel('gauss', fun=fun).get_smoothing(x) - if (y==True).any(): - hs2 = Kernel('gauss', fun=fun).get_smoothing(x[y==True]) - else: - hs2 = 4*hs1 - hsmax = sqrt(hs1*hs2) - for hi in np.linspace(hsmax*0.25,hsmax,9): - plt.figure(k) - k +=1 - unused = kreg_demo3(x,y,fun1, hs=hi, fun=fun, plotlog=False) - - #kreg_demo2(n=n,symmetric=True,fun='hste', plotlog=False) - fig.tile(range(k0,k)) - plt.ioff() - plt.show() - - #test_smoothn_2d() - #test_smoothn_cardioid() \ No newline at end of file diff --git a/pywafo/src/wafo/wafodata.py b/pywafo/src/wafo/wafodata.py index 2440903..3977972 100644 --- a/pywafo/src/wafo/wafodata.py +++ b/pywafo/src/wafo/wafodata.py @@ -86,8 +86,7 @@ class WafoData(object): plotbackend.hold('on') tmp = [] child_args = args if len(args) else tuple(self.plot_args_children) - child_kwds = dict() - child_kwds.update(self.plot_kwds_children) + child_kwds = dict(self.plot_kwds_children).copy() child_kwds.update(**kwds) for child in self.children: tmp1 = child.plot(*child_args, **kwds) @@ -96,8 +95,7 @@ class WafoData(object): if len(tmp) == 0: tmp = None main_args = args if len(args) else tuple(self.plot_args) - main_kwds = dict() - main_kwds.update(self.plot_kwds) + main_kwds = dict(self.plot_kwds).copy() main_kwds.update(kwds) tmp2 = self.plotter.plot(self, *main_args, **main_kwds) return tmp2, tmp @@ -203,7 +201,7 @@ class Plotter_1d(object): def _plot(self, plotflag, wdata, *args, **kwds): x = wdata.args data = transformdata(x, wdata.data, plotflag) - dataCI = () + dataCI = getattr(wdata, 'dataCI', ()) h1 = plot1d(x, data, dataCI, plotflag, *args, **kwds) return h1 @@ -219,7 +217,7 @@ def plot1d(args, data, dataCI, plotflag, *varargin, **kwds): elif plottype == 3: H = plotbackend.stem(args, data, *varargin, **kwds) elif plottype == 4: - H = plotbackend.errorbar(args, data, dataCI[:, 0] - data, dataCI[:, 1] - data, *varargin, **kwds) + H = plotbackend.errorbar(args, data, dataCI[:,0] - data, dataCI[:,1] - data, *varargin, **kwds) elif plottype == 5: H = plotbackend.bar(args, data, *varargin, **kwds) elif plottype == 6: @@ -228,6 +226,9 @@ def plot1d(args, data, dataCI, plotflag, *varargin, **kwds): H = plotbackend.fill_between(args, data, level, *varargin, **kwds); else: H = plotbackend.fill_between(args, data, *varargin, **kwds); + elif plottype==7: + H = plotbackend.plot(args, data, *varargin, **kwds) + H = plotbackend.fill_between(args, dataCI[:,0], dataCI[:,1], alpha=0.2, color='r'); scale = plotscale(plotflag) logXscale = 'x' in scale @@ -255,9 +256,8 @@ def plot1d(args, data, dataCI, plotflag, *varargin, **kwds): plotbackend.axis(ax) - if dataCI and plottype < 3: + if np.any(dataCI) and plottype < 3: plotbackend.hold('on') - plot1d(args, dataCI, (), plotflag, 'r--'); return H