danh
/
pywafo
1
0
Fork 0
Code Issues Pull Requests Releases Wiki Activity

Browse Source
master
Per.Andreas.Brodtkorb 13 years ago
parent c69f54dc0c
commit 8c5754e30c
2 changed files with 82 additions and 99 deletions
Show Stats Download Patch File Download Diff File

125
pywafo/src/wafo/kdetools.py
Unescape Escape View File

@ -1450,7 +1450,7 @@ class Kernel(object):
break break
else: else:
ai = bi ai = bi
y = np.asarray([fun(j) for j in x]) #y = np.asarray([fun(j) for j in x])
#pylab.figure(1) #pylab.figure(1)
#pylab.plot(x,y) #pylab.plot(x,y)
#pylab.show() #pylab.show()
@ -1459,9 +1459,6 @@ class Kernel(object):
try: try:
t_star = optimize.brentq(fun, a=ai, b=bi) t_star = optimize.brentq(fun, a=ai, b=bi)
except: except:
# try:
# t_star = optimize.bisect(fun, a=ai, b=bi+1)
# except:
t_star = 0.28*N**(-2./5) t_star = 0.28*N**(-2./5)
warnings.warn('Failure in obtaining smoothing parameter') warnings.warn('Failure in obtaining smoothing parameter')
@ -2635,68 +2632,19 @@ def kde_demo3():
pylab.figure(0) 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 def kde_demo4(N=50):
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):
'''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior '''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) compared to normal reference rules (in this case the hns)
''' '''
import scipy.stats as st 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) #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.plot(x + loc, st.norm.pdf(x, 0, scale=1)/2, 'k:', label='True density')
pylab.legend() pylab.legend()
def kde_demo6(N=500): def kde_demo5(N=500):
'''Demonstrate that the improved Sheather-Jones plug-in (hisj) is superior '''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) compared to normal reference rules (in this case the hns)
''' '''
import scipy.stats as st import scipy.stats as st
@ -2743,10 +2691,63 @@ def kde_demo6(N=500):
f1.plot() f1.plot()
pylab.plot(data[0], data[1], '.') 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(): def test_docstrings():
import doctest import doctest
doctest.testmod() doctest.testmod()
if __name__ == '__main__': if __name__ == '__main__':
#test_docstrings() #test_docstrings()
kde_demo2() #kde_demo2()
kreg_demo1()

52
pywafo/src/wafo/stats/core.py
Unescape Escape View File

@ -735,33 +735,13 @@ class RegLogit(object):
See also regglm, reglm, regnonlm 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
#% <http://www.gnu.org/licenses/>.
#% Original for MATLAB written by Gordon K Smyth <gks@maths.uq.oz.au>, #% Original for MATLAB written by Gordon K Smyth <gks@maths.uq.oz.au>,
#% U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3, #% U of Queensland, Australia, on Nov 19, 1990. Last revision Aug 3,
#% 1992. #% 1992.
# #
#% Author: Gordon K Smyth <gks@maths.uq.oz.au>, #% Author: Gordon K Smyth <gks@maths.uq.oz.au>,
#% Adapted-By: KH <Kurt.Hornik@wu-wien.ac.at>
#% Revised by: pab #% Revised by: pab
#% -renamed from logistic_regression to reglogit #% -renamed from oridinal to reglogit
#% -added predict, summary and compare #% -added predict, summary and compare
#% Description: Ordinal logistic regression #% Description: Ordinal logistic regression
# #
@ -773,7 +753,7 @@ class RegLogit(object):
self.maxiter =maxiter self.maxiter =maxiter
self.accuracy = accuracy self.accuracy = accuracy
self.alpha = alpha self.alpha = alpha
self.deletcolinear = deletecolinear self.deletecolinear = deletecolinear
self.verbose = False self.verbose = False
self.family = None self.family = None
self.link = None self.link = None
@ -825,8 +805,8 @@ class RegLogit(object):
if len(ix): if len(ix):
X = X[:, ix] X = X[:, ix]
txt = [' %d,' % i for i in iy] txt = [' %d,' % i for i in iy]
txt[-1] = ' %d' % iy[-1] #txt[-1] = ' %d' % iy[-1]
warnings.warn('Covariate matrix is singular. Removing column(s):%s',txt) warnings.warn('Covariate matrix is singular. Removing column(s):%s' % txt)
mx = X.shape[0] mx = X.shape[0]
if (mx != my): if (mx != my):
raise ValueError('x and y must have the same number of observations'); raise ValueError('x and y must have the same number of observations');
@ -872,7 +852,7 @@ class RegLogit(object):
[mx, nx] = X.shape [mx, nx] = X.shape
[my, ny] = y.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)) theta00 = np.log(g / (1 - g))
beta00 = np.zeros((nx, 1)) beta00 = np.zeros((nx, 1))
# starting values # starting values
@ -943,7 +923,7 @@ class RegLogit(object):
else: else:
eta = (y * 0 + 1) * theta.T; eta = (y * 0 + 1) * theta.T;
#end #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) k0 = min(y)
mu = (k0-1)+np.dot(gammai,np.arange(1,nz+2).T); #% E(Y|X) mu = (k0-1)+np.dot(gammai,np.arange(1,nz+2).T); #% E(Y|X)
r = np.corrcoef(np.hstack((y,mu))) r = np.corrcoef(np.hstack((y,mu)))
@ -1234,13 +1214,13 @@ class RegLogit(object):
# first derivative # first derivative
v = g * (1 - g) / p; v = g * (1 - g) / p;
v1 = g1 * (1 - g1) / 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 dl = np.sum(dlogp, axis=0).T
# second derivative # second derivative
w = v * (1 - 2 * g) w = v * (1 - 2 * g)
w1 = v1 * (1 - 2 * g1) 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 return dev, p, dl, d2l
#end %function #end %function
@ -1301,21 +1281,23 @@ def _test_reslife():
mrl.plot() mrl.plot()
def test_reglogit(): def test_reglogit():
y=np.array([1, 1, 2, 1, 3, 2, 3, 2, 3, 3]) y=np.array([1, 1, 2, 1, 3, 2, 3, 2, 3, 3]).reshape(-1,1)
x = np.arange(10).T x = np.arange(10).reshape(-1,1)
b = reglogit(y,x) b = RegLogit()
b.display() % members and methods b.fit(y,x)
#b.display() #% members and methods
b.summary() b.summary()
[mu,plo,pup] = b.predict(); [mu,plo,pup] = b.predict();
plot(x,mu,'g',x,plo,'r:',x,pup,'r:') plot(x,mu,'g',x,plo,'r:',x,pup,'r:')
def main():
def test_doctstrings():
#_test_dispersion_idx() #_test_dispersion_idx()
import doctest import doctest
doctest.testmod() doctest.testmod()
if __name__ == '__main__': if __name__ == '__main__':
pass test_reglogit()
main() #test_doctstrings()(

Write Preview
Loading…
Cancel
Save