gmane.comp.python.scientific.devel

Re: sum_angle() and sum_polar() functions

Stéfan van der Walt — 2012-05-23T23:41:42

Looking at this again, I don't think the comparison is accurate.  The
skimage version of the radon transform does linear interpolation,
whereas sum_angle uses the nearest neighbor.

Stéfan

Adding t-test with unequal variances to stats.py

Junkshops — 2012-05-23T07:38:34

Hi all,

I've issued a pull request (http://github.com/scipy/scipy/pull/227) for 
a version of scipy/stats/stats.py with the following changes:

1) Adds a method for running a t-test with unequal or unknown population 
variances. ttest_ind assumes that population variances are equal.
2) Refactored common code in the 4 t-test methods into shared methods.
3) This section of code, which has variations in multiple methods, looks 
buggy to me:

d = np.mean(a,axis) - np.mean(b,axis)
svar = ((n1-1)*v1+(n2-1)*v2) / float(df)

t = d/np.sqrt(svar*(1.0/n1 + 1.0/n2))
t = np.where((d==0)*(svar==0), 1.0, t) #define t=0/0 = 0, identical means

Surely if d=0, regardless of svar, t should be set to 0, not 1. 
Similarly, if svar = 0 then both variances are zero (assuming that each 
data set has at least 2 points - perhaps there should be a check for 
this?). In that case, if d==0 t should be zero. Otherwise, t should be 
+/-inf. Hence, (svar==0) is redundant.

Accordingly, I've changed the lines in all functions to be the equivalent of

t = np.where((d==0), 0.0, t)

This handles the case where both d and svar are 0. The respective tests 
have also been changed.

If I'm missing something here, please let me know.

Thanks, Gavin

Re: scipy.test() runs 0 tests

Junkshops — 2012-05-22T19:44:02

Hi Skipper,


That's the fix, and the reason the tests are executable is because my 
git repo is in a VirtualBox shared ntfs drive, so Linux permissions 
don't work and every single file in the repo is executable.

I think that if I do

|git config core.filemode false

|

I should be in good shape when I push, although if I'm wrong I'd 
definitely appreciate it if someone sets me right.

Thanks for the help,

Gavin

On 5/22/2012 11:20 AM, Skipper Seabold wrote:
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Re: scipy.test() runs 0 tests

Skipper Seabold — 2012-05-22T18:20:09

How did you build scipy? This works fine for me, though you'll want to
run the tests from the command line since you're not putting scipy on
the path.

[~/src/scipy-skipper] (master)
|11 $ nosetests scipy
<snip>
----------------------------------------------------------------------
Ran 6116 tests in 387.441s

FAILED (SKIP=28, errors=16, failures=19)

Otherwise, you might try doing

nosetest --exe scipy

though I can't imagine why the tests would have the executable bit set
in your source directory.

scipy.test() runs 0 tests

Junkshops — 2012-05-22T15:39:02

Hi all,

Hopefully I'm not doing something too stupid here, but I've googled 
around quite a bit and can't seem to find any tips. Perhaps my google-fu 
is not strong enough. Anyway, I built scipy 0.11.0 in place, seemingly 
successfully, in ~/git/scipy. Ubuntu version is 12.04. I haven't added 
scipy 0.11 to the pypath yet so I'm running in the parent dir.

js< at >ubuntuVB12:~/git/scipy$ sudo aptitude install python-nose
<snip>
Setting up python-nose (1.1.2-3) ...

js< at >ubuntuVB12:~/git/scipy$ python
Python 2.7.3 (default, Apr 20 2012, 22:39:59)
[GCC 4.6.3] on linux2
Type "help", "copyright", "credits" or "license" for more information.
 >>> import scipy
 >>> scipy.test()
Running unit tests for scipy
NumPy version 1.6.1
NumPy is installed in /usr/lib/python2.7/dist-packages/numpy
SciPy version 0.11.0.dev-3852ce2
SciPy is installed in scipy
Python version 2.7.3 (default, Apr 20 2012, 22:39:59) [GCC 4.6.3]
nose version 1.1.2

----------------------------------------------------------------------
Ran 0 tests in 0.161s

OK
<nose.result.TextTestResult run=0 errors=0 failures=0>
 >>> scipy.test("full")
<snip>

----------------------------------------------------------------------
Ran 0 tests in 0.139s

OK
<nose.result.TextTestResult run=0 errors=0 failures=0>

Anyone have any idea what I'm doing wrong?

Thanks very much, Gavin

Re: Code question about signal.ltisys.py

Warren Weckesser — 2012-05-22T10:05:20



Ah, right.  Thanks.

Warren




_______________________________________________
SciPy-Dev mailing list
SciPy-Dev< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Re: warnings in scipy.stats.entropy

Pauli Virtanen — 2012-05-22T09:15:11

[clip]

I'd rather add `xlogy(x, y) = x*log(y)` that treats the case `x==0` specially to
scipy.special, than kludge around the issue by masking or turning warnings off.

Re: Code question about signal.ltisys.py

Pauli Virtanen — 2012-05-22T09:11:56

[clip]

These lines are equivalent, *unless* the class has __setattr__ defined.
That's the case here, so rewriting may end up in different behavior
than intended, depending on what __setattr__ does.

    Pauli

Re: warnings in scipy.stats.entropy

2012-05-22T03:29:19

I works for 2d or nd if qk=None, and uses the (sometimes hidden) default axis=0.

If qk is given, it doesn't work but still uses axis=0 in the sum.

I would say typical state for a stats function that hasn't been
cleaned up. For the ones that I did clean up, I usually added the axis
keyword in cases like this.

Josef

Re: warnings in scipy.stats.entropy

Skipper Seabold — 2012-05-22T02:43:00

Hmm, I didn't think it was intended for 2d cases since there is no
axis keyword and no tests for this. Docstring is unclear, but I've
only used it for 1d and...

import numpy as np

p = np.random.random((10,4))
p[2,3] = 0
q = np.random.random((10,4))
q[2,3] = 0

p /= p.sum(0)
q /= q.sum(0)

from scipy import stats

# bad logic for > 1d
# plus it would return inf, not a 1d array
stats.entropy(p,q)

stats.entropy(p.flatten(), q.flatten())

# len check not shape
q = np.random.random((10,3))

stats.entropy(p, q)

Code question about signal.ltisys.py

Warren Weckesser — 2012-05-22T00:33:06

While taking a look at https://github.com/scipy/scipy/pull/225, which makes
a small changes to signal/ltisys.py,I noticed the surrounding code, which
has lines like this:

            self.__dict__['num'], self.__dict__['den'] = normalize(*args)


My inclination is to rewrite that as

            self.num, self.den = normalize(*args)


But am I missing something?


Warren
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Re: warnings in scipy.stats.entropy

2012-05-21T23:34:54

won't work as replacement, if qk is None then the function is
vectorized for axis=0

array([[ 0.13878479,  0.03527334,  0.12000785,  0.14706888],
       [ 0.07682377,  0.12749588,  0.15172758,  0.19499206],
       [ 0.10462715,  0.1766166 ,  0.        ,  0.09346067],
       [ 0.02208519,  0.14443609,  0.11331574,  0.15090141],
       [ 0.00830154,  0.06009464,  0.05424912,  0.11603281],
       [ 0.05205531,  0.0792505 ,  0.02387006,  0.0061777 ],
       [ 0.00526626,  0.08439299,  0.17298407,  0.09992403],
       [ 0.16510456,  0.07008839,  0.01962196,  0.07101189],
       [ 0.23265325,  0.15908956,  0.2072021 ,  0.08105922],
       [ 0.19429818,  0.06326201,  0.13702153,  0.03937134]])

array([ 1.9678332 ,  2.19817097,  2.0136922 ,  2.1379255 ])

8.3176218626994789
8.3176218626994789

Josef

Re: warnings in scipy.stats.entropy

Skipper Seabold — 2012-05-21T23:23:33

https://github.com/scipy/scipy/pull/226

Re: warnings in scipy.stats.entropy

2012-05-21T23:11:14

+1,

buggy: if qk is given, then the function isn't vectorized.

Josef

Re: warnings in scipy.stats.entropy

Nathaniel Smith — 2012-05-21T22:43:13

I like the mask version better.

- N

warnings in scipy.stats.entropy

Skipper Seabold — 2012-05-21T22:39:43

Currently in scipy.stats.entropy if you are not ignoring them you will
see warnings when the function is given a probability of zero even
though the case of zero is specifically handled in the function.
Rightly or wrongly this makes me cringe. What do people think about
fixing this by using seterr explicitly in the function or masking the
zeros. Eg.,

import numpy as np
from scipy.stats import entropy

prob = np.random.uniform(0,20, size=10)
prob[5] = 0
prob = prob/prob.sum()

np.seterr(all = 'warn')
entropy(prob) # too loud

Instead we could do (within entropy)

oldstate = np.geterr()
np.seterr(divide='ignore', invalid='ignore')
entropy(prob)
np.seterr(**oldstate)

or just mask the zeros in the first place if this is too much

idx = prob > 0
-np.sum(prob[idx] * np.log(prob[idx]))

Thoughts?

Skipper

Faster Hough/Radon transform (was sum_angle()sum_polar() functions)

Charles R Harris — 2012-05-21T21:25:23

Thought I'd advertize this under a different subject heading. The sum_angle
function is really the inner loop but the timings are about 10x faster than
the Radon transform in skimage, so some here might be interested.

Chuck
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Re: scipy.signal.normalize problems?

Josh Lawrence — 2012-05-21T19:11:56

If I change the allclose lines to have

allclose(0, outb[:,0], atol=1e-14) 

it works. I think that is what the original goal was, anyways. From the documentation of allclose, what I have written above result in ensuring

np.abs(0 - outb[:,0]) > atol + rtol * np.abs(outb[:,0])

with rtol defaulting to 1e-5. I'm still not sure about how things have been written for the 'b' argument of normalize being rank-2, so I can't guarantee that my fix makes things work for that. Should I file a bug report/submit a patch/send a git pull request, etc?

Cheers,

Josh Lawrence

On May 21, 2012, at 2:45 PM, Josh Lawrence wrote:

scipy.signal.normalize problems?

Josh Lawrence — 2012-05-21T18:45:55

Hey all,

I've been having some problems designing a Chebyshev filter and I think I have narrowed down the hang-up to scipy.signal.normalize. I think what's going on in my case is that line 286 of filter_design.py (the first allclose call in the normalize function) is producing a false positive. Here's the function definition:

def normalize(b, a):
    """Normalize polynomial representation of a transfer function.

    If values of b are too close to 0, they are removed. In that case, a
    BadCoefficients warning is emitted.
    """
    b, a = map(atleast_1d, (b, a))
    if len(a.shape) != 1:
        raise ValueError("Denominator polynomial must be rank-1 array.")
    if len(b.shape) > 2:
        raise ValueError("Numerator polynomial must be rank-1 or"
                         " rank-2 array.")
    if len(b.shape) == 1:
        b = asarray([b], b.dtype.char)
    while a[0] == 0.0 and len(a) > 1:
        a = a[1:]
    outb = b * (1.0) / a[0]
    outa = a * (1.0) / a[0]
    if allclose(outb[:, 0], 0, rtol=1e-14): <------------------ Line 286
        warnings.warn("Badly conditioned filter coefficients (numerator): the "
                      "results may be meaningless", BadCoefficients)
        while allclose(outb[:, 0], 0, rtol=1e-14) and (outb.shape[-1] > 1):
            outb = outb[:, 1:]
    if outb.shape[0] == 1:
        outb = outb[0]
    return outb, outa

I marked line 286. If I reproduce all the steps carried out by scipy.signal.iirdesign, I end up with a (b, a) pair which results of scipy.signal.lp2lp and looks like this:

In [106]: b_lp2
Out[106]: array([  1.55431359e-06+0.j])

In [107]: a_lp2
Out[107]: 
array([  1.00000000e+00 +0.00000000e+00j,
         3.46306104e-01 -2.01282794e-16j,
         2.42572185e-01 -6.08207573e-17j,
         5.92946943e-02 +0.00000000e+00j,
         1.82069156e-02 +5.55318531e-18j,
         2.89328123e-03 +0.00000000e+00j,
         4.36566281e-04 -2.95766719e-19j,
         3.50842810e-05 -3.19180568e-20j,   1.64641246e-06 -1.00966301e-21j])

scipy.signal.iirdesign takes b_lp2, a_lp2 (my local variable names to keep track of what's going on) and runs them through scipy.signal.bilinear (in filter_design.py bilinear is called on line 624 within iirfilter. iirdesign calls iirfilter which calls bilinear). Inside bilinear, normalize is called on line 445. I've made my own class with bilinear copied and pasted from filter_design.py to test things. In bilinear, the input to normalize is given by

b = [  1.55431359e-06   1.24345087e-05   4.35207804e-05   8.70415608e-05
   1.08801951e-04   8.70415608e-05   4.35207804e-05   1.24345087e-05
   1.55431359e-06]
a = [   72269.02590913  -562426.61430468  1918276.19173089 -3745112.83646825
  4577612.13937628 -3586970.61385926  1759651.18184723  -494097.93515708
    60799.46134722]

In normalize, right before the allclose() call, outb is defined by

outb = [[  2.15073272e-11   1.72058618e-10   6.02205162e-10   1.20441032e-09
    1.50551290e-09   1.20441032e-09   6.02205162e-10   1.72058618e-10
    2.15073272e-11]]


b =

   1.0e-08 *

  Columns 1 through 5

   0.002150733144728   0.017205865157826   0.060220528052392   0.120441056104784   0.150551320130980

  Columns 6 through 9

   0.120441056104784   0.060220528052392   0.017205865157826   0.002150733144728

which matches up rather well for several significant figures.

I apologize if this is not clearly explained, but I'm not sure what to do. I tried messing around with the arguments to allclose (switching it to be allclose(0, outb[:,0], ...), or changing the keyword from rtol to atol). I am also not sure why normalize is setup to run on rank-2 arrays. I looked through filter_design.py and none of the functions contained in it send a rank-2 array to normalize from what I can tell. Any thoughts?

Cheers,

Josh Lawrence

Re: sum_angle() and sum_polar() functions

Charles R Harris — 2012-05-21T16:20:40

Looks like you need to sell this as a faster Radon or Hough, angle_sum is a
bit too generic ;)

Chuck
_______________________________________________
SciPy-Dev mailing list
SciPy-Dev< at >scipy.org
http://mail.scipy.org/mailman/listinfo/scipy-dev

Re: sum_angle() and sum_polar() functions

Robert Jördens — 2012-05-21T08:06:29

On Mon, May 21, 2012 at 1:16 AM, Ralf Gommers
<ralf.gommers< at >googlemail.com> wrote:

In [47]: from scipy.misc.pilutil import radon as scipy_radon
In [50]: from skimage.transform import radon as scikits_radon

In [46]: a=np.random.randn(1000, 1000)

In [53]: angle_sum(a, .1).shape
Out[53]: (1100,)

In [54]: scikits_radon(a, [.1]).shape
Out[54]: (1415, 1)

In [55]: scipy_radon(a, [.1]).shape
Out[55]: (1000, 1)

In [56]: %timeit angle_sum(a, .1)
10 loops, best of 3: 24.7 ms per loop

In [57]: %timeit scikits_radon(a, [.1])
10 loops, best of 3: 136 ms per loop

In [58]: %timeit scipy_radon(a, [.1])
10 loops, best of 3: 47.6 ms per loop

In [59]: a.sum()
Out[59]: -682.83728184031156

In [60]: angle_sum(a, .1).sum()
Out[60]: -682.83728184032327

In [61]: scikits_radon(a, [.1]).sum()
Out[61]: -682.83727334014043

In [62]: scipy_radon(a, [.1]).sum()
Out[62]: 124316062.0


angle_sum() at least can be applied beyond image processing in the
analysis of coupling matrices or adjacency graphs. Both functions are
simple and generic.