From 8c5754e30c88802debbd5954e300fe45f016bf33 Mon Sep 17 00:00:00 2001 From: "Per.Andreas.Brodtkorb" Date: Mon, 28 Nov 2011 15:44:44 +0000 Subject: [PATCH] --- pywafo/src/wafo/kdetools.py | 127 +++++++++++++++++----------------- pywafo/src/wafo/stats/core.py | 54 +++++---------- 2 files changed, 82 insertions(+), 99 deletions(-) diff --git a/pywafo/src/wafo/kdetools.py b/pywafo/src/wafo/kdetools.py index 09c5de6..6a5a864 100644 --- a/pywafo/src/wafo/kdetools.py +++ b/pywafo/src/wafo/kdetools.py @@ -1450,7 +1450,7 @@ class Kernel(object): break else: ai = bi - y = np.asarray([fun(j) for j in x]) + #y = np.asarray([fun(j) for j in x]) #pylab.figure(1) #pylab.plot(x,y) #pylab.show() @@ -1459,9 +1459,6 @@ class Kernel(object): try: t_star = optimize.brentq(fun, a=ai, b=bi) except: -# try: -# t_star = optimize.bisect(fun, a=ai, b=bi+1) -# except: t_star = 0.28*N**(-2./5) warnings.warn('Failure in obtaining smoothing parameter') @@ -2635,68 +2632,19 @@ def kde_demo3(): pylab.figure(0) -def kde_demo4(hs=None, fast=False): - ''' - - ''' - N = 100 - #ei = np.random.normal(loc=0, scale=0.075, size=(N,)) - ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525, - 0.06096313, -0.16572389, 0.01838912, -0.06251845, -0.09186784, - -0.04304887, -0.13365788, -0.0185279 , -0.07289167, 0.02319097, - 0.06887854, -0.08938374, -0.15181813, 0.03307712, 0.08523183, - -0.0378058 , -0.06312874, 0.01485772, 0.06307944, -0.0632959 , - 0.18963205, 0.0369126 , -0.01485447, 0.04037722, 0.0085057 , - -0.06912903, 0.02073998, 0.1174351 , 0.17599277, -0.06842139, - 0.12587608, 0.07698113, -0.0032394 , -0.12045792, -0.03132877, - 0.05047314, 0.02013453, 0.04080741, 0.00158392, 0.10237899, - -0.09069682, 0.09242174, -0.15445323, 0.09190278, 0.07138498, - 0.03002497, 0.02495252, 0.01286942, 0.06449978, 0.03031802, - 0.11754861, -0.02322272, 0.00455867, -0.02132251, 0.09119446, - -0.03210086, -0.06509545, 0.07306443, 0.04330647, 0.078111 , - -0.04146907, 0.05705476, 0.02492201, -0.03200572, -0.02859788, - -0.05893749, 0.00089538, 0.0432551 , 0.04001474, 0.04888828, - -0.17708392, 0.16478644, 0.1171006 , 0.11664846, 0.01410477, - -0.12458953, -0.11692081, 0.0413047 , -0.09292439, -0.07042327, - 0.14119701, -0.05114335, 0.04994696, -0.09520663, 0.04829406, - -0.01603065, -0.1933216 , 0.19352763, 0.11819496, 0.04567619, - -0.08348306, 0.00812816, -0.00908206, 0.14528945, 0.02901065]) - x = np.linspace(0, 1, N) - - - y0 = 2*np.exp(-x**2/(2*0.3**2))+3*np.exp(-(x-1)**2/(2*0.7**2)) - y = y0 + ei - kreg = KRegression(x, y, p=0, hs=hs) - kreg.tkde.kernel.get_smooting = kreg.tkde.kernel.hste - if fast: - kreg.__call__ = kreg.eval_grid_fast - - f = kreg(output='plot', title='Kernel regression', plotflag=1) - pylab.figure(0) - f.plot(label='p=0') - - - kreg.p=1 - f1 = kreg(output='plot', title='Kernel regression', plotflag=1) - f1.plot(label='p=1') - pylab.plot(x,y,'.', x,y0, 'k') - pylab.legend() - - pylab.show() - - print(kreg.tkde.tkde.inv_hs) - print(kreg.tkde.tkde.hs) -def kde_demo5(N=50): + +def kde_demo4(N=50): '''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior - for multimodal distributions + for 1D multimodal distributions - KDEDEMO5 shows that the improved Sheather-Jones plug-in smoothing is a better + KDEDEMO4 shows that the improved Sheather-Jones plug-in smoothing is a better compared to normal reference rules (in this case the hns) ''' import scipy.stats as st - data = np.hstack((st.norm.rvs(loc=5, scale=1, size=(N,)), st.norm.rvs(loc=-5, scale=1, size=(N,)))) + data = np.hstack((st.norm.rvs(loc=5, scale=1, size=(N,)), + st.norm.rvs(loc=-5, scale=1, size=(N,)))) #x = np.linspace(1.5e-3, 5, 55) @@ -2717,11 +2665,11 @@ def kde_demo5(N=50): pylab.plot(x + loc, st.norm.pdf(x, 0, scale=1)/2, 'k:', label='True density') pylab.legend() -def kde_demo6(N=500): +def kde_demo5(N=500): '''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior - for multimodal distributions + for 2D multimodal distributions - KDEDEMO5 shows that the improved Sheather-Jones plug-in smoothing is a better + KDEDEMO5 shows that the improved Sheather-Jones plug-in smoothing is better compared to normal reference rules (in this case the hns) ''' import scipy.stats as st @@ -2742,11 +2690,64 @@ def kde_demo6(N=500): pylab.clf() f1.plot() pylab.plot(data[0], data[1], '.') + +def kreg_demo1(hs=None, fast=False, fun='hisj'): + ''' + + ''' + N = 100 + #ei = np.random.normal(loc=0, scale=0.075, size=(N,)) + ei = np.array([-0.08508516, 0.10462496, 0.07694448, -0.03080661, 0.05777525, + 0.06096313, -0.16572389, 0.01838912, -0.06251845, -0.09186784, + -0.04304887, -0.13365788, -0.0185279 , -0.07289167, 0.02319097, + 0.06887854, -0.08938374, -0.15181813, 0.03307712, 0.08523183, + -0.0378058 , -0.06312874, 0.01485772, 0.06307944, -0.0632959 , + 0.18963205, 0.0369126 , -0.01485447, 0.04037722, 0.0085057 , + -0.06912903, 0.02073998, 0.1174351 , 0.17599277, -0.06842139, + 0.12587608, 0.07698113, -0.0032394 , -0.12045792, -0.03132877, + 0.05047314, 0.02013453, 0.04080741, 0.00158392, 0.10237899, + -0.09069682, 0.09242174, -0.15445323, 0.09190278, 0.07138498, + 0.03002497, 0.02495252, 0.01286942, 0.06449978, 0.03031802, + 0.11754861, -0.02322272, 0.00455867, -0.02132251, 0.09119446, + -0.03210086, -0.06509545, 0.07306443, 0.04330647, 0.078111 , + -0.04146907, 0.05705476, 0.02492201, -0.03200572, -0.02859788, + -0.05893749, 0.00089538, 0.0432551 , 0.04001474, 0.04888828, + -0.17708392, 0.16478644, 0.1171006 , 0.11664846, 0.01410477, + -0.12458953, -0.11692081, 0.0413047 , -0.09292439, -0.07042327, + 0.14119701, -0.05114335, 0.04994696, -0.09520663, 0.04829406, + -0.01603065, -0.1933216 , 0.19352763, 0.11819496, 0.04567619, + -0.08348306, 0.00812816, -0.00908206, 0.14528945, 0.02901065]) + x = np.linspace(0, 1, N) + + + y0 = 2*np.exp(-x**2/(2*0.3**2))+3*np.exp(-(x-1)**2/(2*0.7**2)) + y = y0 + ei + kernel = Kernel('gauss',fun=fun) + kreg = KRegression(x, y, p=0, hs=hs, kernel=kernel) + if fast: + kreg.__call__ = kreg.eval_grid_fast + + f = kreg(output='plot', title='Kernel regression', plotflag=1) + pylab.figure(0) + f.plot(label='p=0') + + + kreg.p=1 + f1 = kreg(output='plot', title='Kernel regression', plotflag=1) + f1.plot(label='p=1') + pylab.plot(x,y,'.', x,y0, 'k') + pylab.legend() + + pylab.show() + print(kreg.tkde.tkde.inv_hs) + print(kreg.tkde.tkde.hs) + def test_docstrings(): import doctest doctest.testmod() if __name__ == '__main__': #test_docstrings() - kde_demo2() \ No newline at end of file + #kde_demo2() + kreg_demo1() \ No newline at end of file diff --git a/pywafo/src/wafo/stats/core.py b/pywafo/src/wafo/stats/core.py index 212a187..4523df7 100644 --- a/pywafo/src/wafo/stats/core.py +++ b/pywafo/src/wafo/stats/core.py @@ -734,34 +734,14 @@ class RegLogit(object): See also regglm, reglm, regnonlm ''' - - #% Copyright (C) 1995, 1996, 1997, 1998, 1999, 2000, 2002, 2005, 2007 - #% Kurt Hornik - #% - #% Reglogit is free software; you can redistribute it and/or modify it - #% under the terms of the GNU General Public License as published by - #% the Free Software Foundation; either version 3 of the License, or (at - #% your option) any later version. - #% - #% Reglogit is distributed in the hope that it will be useful, but - #% WITHOUT ANY WARRANTY; without even the implied warranty of - #% MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU - #% General Public License for more details. - #% - #% You should have received a copy of the GNU General Public License - #% along with Reglogit; see the file COPYING. If not, see - #% . - - - + #% Original for MATLAB written by Gordon K Smyth , #% U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, #% 1992. # #% Author: Gordon K Smyth , - #% Adapted-By: KH #% Revised by: pab - #% -renamed from logistic_regression to reglogit + #% -renamed from oridinal to reglogit #% -added predict, summary and compare #% Description: Ordinal logistic regression # @@ -773,7 +753,7 @@ class RegLogit(object): self.maxiter =maxiter self.accuracy = accuracy self.alpha = alpha - self.deletcolinear = deletecolinear + self.deletecolinear = deletecolinear self.verbose = False self.family = None self.link = None @@ -825,8 +805,8 @@ class RegLogit(object): if len(ix): X = X[:, ix] txt = [' %d,' % i for i in iy] - txt[-1] = ' %d' % iy[-1] - warnings.warn('Covariate matrix is singular. Removing column(s):%s',txt) + #txt[-1] = ' %d' % iy[-1] + warnings.warn('Covariate matrix is singular. Removing column(s):%s' % txt) mx = X.shape[0] if (mx != my): raise ValueError('x and y must have the same number of observations'); @@ -872,7 +852,7 @@ class RegLogit(object): [mx, nx] = X.shape [my, ny] = y.shape - g = z.sum(axis=0).cumsum() / my + g = (z.sum(axis=0).cumsum() / my).reshape(-1,1) theta00 = np.log(g / (1 - g)) beta00 = np.zeros((nx, 1)) # starting values @@ -943,7 +923,7 @@ class RegLogit(object): else: eta = (y * 0 + 1) * theta.T; #end - gammai = np.diff(np.hstack(((y * 0), self.logitinv(eta), (y * 0 + 1))),1,2) + gammai = np.diff(np.hstack(((y * 0), self.logitinv(eta), (y * 0 + 1))),n=1,axis=1) k0 = min(y) mu = (k0-1)+np.dot(gammai,np.arange(1,nz+2).T); #% E(Y|X) r = np.corrcoef(np.hstack((y,mu))) @@ -1234,13 +1214,13 @@ class RegLogit(object): # first derivative v = g * (1 - g) / p; v1 = g1 * (1 - g1) / p; - dlogp = np.hstack(((dmult(v, z) - dmult(v1, z1)), (dmult(v - v1, x)))) + dlogp = np.hstack((((v*z) - (v1*z1)), ((v - v1)*x))) dl = np.sum(dlogp, axis=0).T # second derivative w = v * (1 - 2 * g) w1 = v1 * (1 - 2 * g1) - d2l = zx.T * dmult (w, zx) - z1x.T * dmult(w1, z1x) - dlogp.T * dlogp; + d2l = np.dot(zx.T, (w*zx)) - np.dot(z1x.T, (w1*z1x)) - np.dot(dlogp.T, dlogp) return dev, p, dl, d2l #end %function @@ -1301,21 +1281,23 @@ def _test_reslife(): mrl.plot() def test_reglogit(): - y=np.array([1, 1, 2, 1, 3, 2, 3, 2, 3, 3]) - x = np.arange(10).T - b = reglogit(y,x) - b.display() % members and methods + y=np.array([1, 1, 2, 1, 3, 2, 3, 2, 3, 3]).reshape(-1,1) + x = np.arange(10).reshape(-1,1) + b = RegLogit() + b.fit(y,x) + #b.display() #% members and methods b.summary() [mu,plo,pup] = b.predict(); plot(x,mu,'g',x,plo,'r:',x,pup,'r:') -def main(): + +def test_doctstrings(): #_test_dispersion_idx() import doctest doctest.testmod() if __name__ == '__main__': - pass - main() + test_reglogit() + #test_doctstrings()(