From 25b47ebc048b050680a77a2e958f47939f73e0a1 Mon Sep 17 00:00:00 2001 From: "per.andreas.brodtkorb" Date: Tue, 27 Dec 2011 22:55:31 +0000 Subject: [PATCH] --- pywafo/src/wafo/kdetools.py | 37 ++++++++++++++++++++++--------------- 1 file changed, 22 insertions(+), 15 deletions(-) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index c1aae43..389f3be 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -3471,12 +3471,12 @@ def _logitinv(x): return 1.0/(np.exp(-x)+1) -def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): +def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj', plotlog=False): import scipy.stats as st from sg_filter import SavitzkyGolay dist = st.norm - scale1 = 0.7 - loc1= 2 + scale1 = 0.3 + loc1= 0 norm1 = dist.pdf(-loc1, loc=-loc1, scale=scale1) + dist.pdf(-loc1, loc=loc1, scale=scale1) fun1 = lambda x : (dist.pdf(x, loc=-loc1, scale=scale1) + dist.pdf(x, loc=loc1, scale=scale1))/norm1 @@ -3501,18 +3501,25 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): y = yi[i] xmin, xmax = x.min(), x.max() - ni = int(2*(xmax-xmin)/hopt) + ni = 2*int((xmax-xmin)/hopt)+1 print(ni) print(xmin, xmax) xi = np.linspace(xmin-hopt,xmax+hopt, ni) - xiii = np.linspace(xmin-hopt,xmax+hopt, 4*ni) + xiii = np.linspace(xmin-hopt,xmax+hopt, 4*ni+1) + + kreg = KRegression(x, y, hs=hs*0.5) + fi = kreg(xi) + f = kreg(xiii,output='plotobj', title='Kernel regression n=%d, %s=%g' % (n,fun,hs), plotflag=1) + c = gridcount(x, xi) c0 = gridcount(x[y==True],xi) yi = np.where(c==0, 0, c0/c) - yi[yi==0] = 1.0/(c.min()+4)#yi[yi>0].min()/sqrt(n) - yi[yi==1] = 1-1.0/(c.min()+4) #(1-yi[yi<1].max())/sqrt(n) + #yi[yi==0] = 1.0/(c[c!=0].min()+4)# + #yi[yi==1] = 1-1.0/(c[c!=0].min()+4) # + yi[yi==0] = fi[yi==0] #yi[yi>0].min()/sqrt(n) + yi[yi==1] =1-(1-yi[yi<1].max())/sqrt(n) logity =_logit(yi) # plt.plot(xi, np.log(yi/(1-yi)), xi,logity,'.') @@ -3521,7 +3528,7 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): logity[logity==-40]=np.nan slogity = smoothn(logity, robust=True) - slogity2 = SavitzkyGolay(n=3, degree=2).smooth(logity) + slogity2 = SavitzkyGolay(n=2, degree=2).smooth(logity) sa1 = sqrt(evar(logity)) sa = (slogity-logity).std() print('estd = %g %g' % (sa,sa1)) @@ -3539,7 +3546,7 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): #return #print(logyi) - gkreg = KRegression(xi, logity, hs=hs/2, xmin=xmin-2*hopt,xmax=xmax+2*hopt) + gkreg = KRegression(xi, logity, hs=hs*0.4, xmin=xmin-2*hopt,xmax=xmax+2*hopt) fg = gkreg.eval_grid(xi,output='plotobj', title='Kernel regression', plotflag=1) sa = (fg.data-logity).std() sa2 = iqrange(fg.data-logity) / 1.349 @@ -3554,7 +3561,6 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): # return - fg = gkreg.eval_grid(xiii,output='plotobj', title='Kernel regression', plotflag=1) dx = xi[1]-xi[0] @@ -3569,6 +3575,9 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): pup3 = _logitinv(fg.data+z0*sa) fg.data = pi + + #pi = f.data + # ref Casella and Berger (1990) "Statistical inference" pp444 a = 2*pi + z0**2/(ciii+1e-16) @@ -3590,9 +3599,6 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): plt.figure(2) fg.plot(label='KReg grid') - - kreg = KRegression(x, y, hs=hs*0.5) - f = kreg(xiii,output='plotobj', title='Kernel regression n=%d, %s=%g' % (n,fun,hs), plotflag=1) f.plot(label='KRegression') labtxt = '%d CI' % (int(100*(1-alpha))) plt.plot(xi, syi, 'k',xi, syi2,'k--', label='smoothn') @@ -3605,7 +3611,8 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'): print(kreg.tkde.tkde.hs) plt.legend() h = plt.gca() - #plt.setp(h,yscale='log') + if plotlog: + plt.setp(h,yscale='log') plt.show() def kde_gauss_demo(n=50): @@ -3672,6 +3679,6 @@ if __name__ == '__main__': #kde_demo2() #kreg_demo1(fast=True) #kde_gauss_demo() - kreg_demo2(n=3800,symmetric=True,fun='hste') + kreg_demo2(n=180,symmetric=True,fun='hste', plotlog=False) #test_smoothn_2d() #test_smoothn_cardioid() \ No newline at end of file