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()(