From 02f8fbf9492dea02d2bead610432d89ad5719163 Mon Sep 17 00:00:00 2001
From: "per.andreas.brodtkorb"
 <per.andreas.brodtkorb@96a810ac-173e-11df-b9d6-e102a254d9bb>
Date: Tue, 27 Dec 2011 21:39:21 +0000
Subject: [PATCH] updated kreg_demo2

---
 pywafo/src/wafo/kdetools.py | 28 +++++++++++++++++++---------
 1 file changed, 19 insertions(+), 9 deletions(-)

diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py
index de32b38..c1aae43 100644
--- a/pywafo/src/wafo/kdetools.py
+++ b/pywafo/src/wafo/kdetools.py
@@ -3475,8 +3475,8 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'):
     import scipy.stats as st
     from sg_filter import SavitzkyGolay
     dist = st.norm
-    scale1 = 0.4
-    loc1= 1
+    scale1 = 0.7
+    loc1= 2
     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 
 
@@ -3511,8 +3511,8 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'):
     c0 = gridcount(x[y==True],xi)
     
     yi = np.where(c==0, 0, c0/c)
-    yi[yi==0] = yi[yi>0].min()/sqrt(n)
-    yi[yi==1] = 1-1.0/n #(1-yi[yi<1].max())/sqrt(n)
+    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)
     
     logity =_logit(yi)
 #    plt.plot(xi, np.log(yi/(1-yi)), xi,logity,'.')
@@ -3539,7 +3539,7 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'):
     #return
     #print(logyi)
     
-    gkreg = KRegression(xi, logity, hs=hs/3.5, xmin=xmin-2*hopt,xmax=xmax+2*hopt)
+    gkreg = KRegression(xi, logity, hs=hs/2, 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
@@ -3577,25 +3577,35 @@ def kreg_demo2(n=100, hs=None, symmetric=False, fun='hisj'):
     pup2 = ((a+sqrt(a**2-2*pi**2*b))/b).clip(min=0,max=1)
     pup = (pi + z0*np.sqrt(pi*(1-pi)/ciii)).clip(min=0,max=1)
     plo = (pi - z0*np.sqrt(pi*(1-pi)/ciii)).clip(min=0,max=1)
+    
+    den = 1+(z0**2./ciii);
+    xc=(pi+(z0**2)/(2*ciii))/den;
+    halfwidth=(z0*sqrt((pi*(1-pi)/ciii)+(z0**2/(4*(ciii**2)))))/den
+    plo = xc-halfwidth # wilson score
+    pup = xc+halfwidth # wilson score
+
+    
     #print(fg.data)
     #fg.data = np.exp(fg.data)
     
     plt.figure(2)
     fg.plot(label='KReg grid')
     
-    kreg = KRegression(x, y, hs=hs)
+    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')
-    plt.fill_between(xiii, pup, plo, alpha=0.15,color='r', linestyle='--', label=labtxt)
-    plt.fill_between(xiii, pup2, plo2,alpha = 0.10, color='b', linestyle=':',label='%d CI2' % (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))))
     #plt.plot(xiii, 0.5*np.cos(xiii)+.5, 'r', label='True model')
     plt.plot(xiii, fun1(xiii), 'r', label='True model')
     plt.scatter(xi,yi, label='data')
     print(np.nanmax(f.data))
     print(kreg.tkde.tkde.hs)
     plt.legend()
+    h = plt.gca()
+    #plt.setp(h,yscale='log')
     plt.show()
     
 def kde_gauss_demo(n=50):
@@ -3662,6 +3672,6 @@ if __name__ == '__main__':
     #kde_demo2()
     #kreg_demo1(fast=True)
     #kde_gauss_demo()
-    kreg_demo2(n=750,symmetric=True,fun='hisj')
+    kreg_demo2(n=3800,symmetric=True,fun='hste')
     #test_smoothn_2d()
     #test_smoothn_cardioid()
\ No newline at end of file