Mailman 3 March 2009 - SciPy-Dev

Possible Error in Kendall's Tau (scipy.stats.stats.kendalltau)
by Almer S. Tigelaar 19 Mar '09

19 Mar '09

Hello, (I realize this mail is a bit lengthy, but I would appreciate it if someone could comment on it). I believe that I found a bug in your implementation of Kendall's Tau. I have evaluated the implementation (to verify a self-written implementation). When the results turned out to be different I investigated the current SciPy implementation at the following URL: http://svn.scipy.org/svn/scipy/trunk/scipy/stats/stats.py (I am aware of the fact that there is also a Kendall's Tau implementation in mstats.py, but have not evaluated that implementation yet). I will give some explanation of my interpretation of Kendall's tau, an example showing the differences between SciPy's and my implementation and a possible fix for SciPy's implementation. Your implementation is Kendall's tau-b with tie correction (same as mine). I take as my reference definition, the one in the following 'poster' paper: http://portal.acm.org/citation.cfm?id=1277935 (this same definition appears in other places as well, this is the shortest resource I could find) Recall that Kendall's tau calculates a score t given two rankings R1 and R2. Variables P, Q and T are all characteristics of the pairs in those rankings. The definition given in the reference is: t = (P - Q) / SQRT((P + Q + T) * (P + Q + U)) where P is the number of concordant pairs, Q the number of discordant pairs, T the number ties in R1 and U the number of ties in R2. An example: ----------- Let's use two identical rankings with a tie: A B C R1 = [1, 1, 2] R2 = [1, 1, 2] There are three pair combinations in these lists, namely: (A, B), (A, C) and (B, C). It is obvious that _one_ of these combinations has a tie for both lists (the (A,B) combination which is (1,1) for both R1 and R2). So, since there is one tie in both list we have T = U = 1 We find that there are two concordant pairs in both lists (A, C) and (B,C) so P = 2. There are no discordant pairs, so Q = 0. With all variables given, we can now calculate Kendall's tau for R1 and R2: t = (2 - 0) / SQRT((2 + 0 + 1)*(2 + 0 + 1)) t = 2 / SQRT(3*3) t = 2 / 3 t = 0.6666666 However, using scipy (svn HEAD) as follows: import scipy.stats.stats as s s.kendalltau([1,1,2], [1,1,2]) Yields t = 1.0: (1.0, 0.11718509694604401) Which I believe is wrong (or at least: has no correction for ties, as is claimed in the source code). If there are three combinations and one of these is a tie, and the other two combinations are concordant, it makes sense that Kendall's tau-b should yield 2 / 3. The cause and fix ----------------- Playing around with SciPy's code (and comparing it with my own) I believe I discovered a probable cause for this difference in SciPy's code. Again, I used the implementation at the following URL: http://svn.scipy.org/svn/scipy/trunk/scipy/stats/stats.py (please take look at the implementation first, otherwise you will not understand my explanation) In the 'kendalltau(x,y)' function we see a test for ties and an 'else' branch. In the 'else' branch the values of 'n1' and 'n2' are incremented if there is a tie (conforming to +T and +U in the formula given above). However, I believe that the 'if' conditions here are wrong: 1) Consider that if 'a1' has value '0' it is tied (the same goes for 'a2'). In the else branch I see: if a1: n1 = n1 + 1 if a2: n2 = n2 + 1 So, here the addition takes places on the variables (n1, n2) if there is NO tie, instead of if there is a tie. Hence, this explains the different outcome. Translating this back to the formula gives me T = U = 0, which would yield: t = (2 - 0) / SQRT((2 + 0 + 0)*(2 + 0 + 0)) t = 2 / SQRT(2*2) t = 2 / 2 t = 1.0 Which is indeed consistent with the SciPy outcome. Henceforth, I believe the solution to this is to correct the condition in the if statements in the Kendall's tau function: if not a1: n1 = n1 + 1 if not a2: n2 = n2 + 1 Closing ------- Of course, my interpretation of Kendall's Tau could be wrong. Since I can not exclude that possibility I would appreciate it if one of you could check and see if you reach the same conclusion. Maybe the base formula that SciPy uses is different. I have compared your implementation also to that implemented in the R project, however their source code suggests that they do not adjust for ties (effectively implementing Kendall's tau-a). -- With kind regards, Almer S. Tigelaar University of Twente

4 34

Status on ubuntu jaunty 64bits.
by Xavier Gnata 17 Mar '09

17 Mar '09

Hi, Here it is: ====================================================================== ERROR: test_implicit (test_odr.TestODR) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/odr/tests/test_odr.py", line 88, in test_implicit out = implicit_odr.run() File "/usr/local/lib/python2.6/dist-packages/scipy/odr/odrpack.py", line 1055, in run self.output = Output(apply(odr, args, kwds)) TypeError: y must be a sequence or integer (if model is implicit) ====================================================================== FAIL: test_random_real (test_basic.TestSingleIFFT) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/fftpack/tests/test_basic.py", line 206, in test_random_real assert_array_almost_equal (y2, x) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 311, in assert_array_almost_equal header='Arrays are not almost equal') File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 296, in assert_array_compare raise AssertionError(msg) AssertionError: Arrays are not almost equal (mismatch 0.900900900901%) x: array([ 0.89560729 -4.65661287e-09j, 0.87991965 +7.21774995e-09j, 0.44631395 -2.04890966e-08j, 0.71974921 +4.46598869e-11j, 0.20776373 +1.41855736e-08j, 0.83089650 -1.69798398e-09j,... y: array([ 0.89560729, 0.87991947, 0.44631401, 0.71974903, 0.20776364, 0.83089662, 0.86079419, 0.93193549, 0.20852582, 0.51215041, 0.91066802, 0.99397069, 0.74227983, 0.67712617, 0.244197 ,... ====================================================================== FAIL: test_iv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1653, in test_iv_cephes_vs_amos self.check_cephes_vs_amos(iv, iv, rtol=1e-8, atol=1e-305) File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1642, in check_cephes_vs_amos assert_tol_equal(c1, c2, err_msg=(v, z), rtol=rtol, atol=atol) File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 296, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-08, atol=1e-305 (-120, -11) (mismatch 100.0%) x: array(1.3384173609003782e-110) y: array((1.3384173859242368e-110+0j)) ====================================================================== FAIL: test_yn_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1598, in test_yn_zeros 488.98055964441374646], rtol=1e-19) File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 296, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-19, atol=0 (mismatch 100.0%) x: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) y: array([ 450.136, 463.057, 472.807, 481.274, 488.981]) ====================================================================== FAIL: test_ynp_zeros (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1604, in test_ynp_zeros assert_tol_equal(yvp(443, ao), 0, atol=1e-15) File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 38, in assert_tol_equal verbose=verbose, header=header) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 296, in assert_array_compare raise AssertionError(msg) AssertionError: Not equal to tolerance rtol=1e-07, atol=1e-15 (mismatch 100.0%) x: array([ 1.239e-10, -8.119e-16, 3.608e-16, 5.898e-16, 1.226e-15]) y: array(0) ====================================================================== FAIL: test_yv_cephes_vs_amos (test_basic.TestBessel) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1650, in test_yv_cephes_vs_amos self.check_cephes_vs_amos(yv, yn, rtol=1e-11, atol=1e-305) File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 1640, in check_cephes_vs_amos assert c2.imag != 0, (v, z) AssertionError: (301, 1.0) ====================================================================== FAIL: test_pbdv (test_basic.TestCephes) ---------------------------------------------------------------------- Traceback (most recent call last): File "/usr/local/lib/python2.6/dist-packages/scipy/special/tests/test_basic.py", line 370, in test_pbdv assert_equal(cephes.pbdv(1,0),(0.0,0.0)) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 176, in assert_equal assert_equal(actual[k], desired[k], 'item=%r\n%s' % (k,err_msg), verbose) File "/usr/lib/python2.6/dist-packages/numpy/testing/utils.py", line 183, in assert_equal raise AssertionError(msg) AssertionError: Items are not equal: item=1 ACTUAL: 1.0 DESIRED: 0.0 ---------------------------------------------------------------------- The last error on pbdv is a quite old one. Can someone reproduce this error? http://projects.scipy.org/scipy/ticket/803 According to mathematica (and my understanding of pbdv....) (0.0, 1.0) is the correct answer. Xavier

2 2

Warning about remaining issues in stats.distributions ?
by josef.pktd＠gmail.com 16 Mar '09

16 Mar '09

There are still some known failures or methods, that don't work very well for some distributions: * skew. kurtosis: a few remaining incorrect results, no exact test available * entropy: returns nan in some cases, correctness of values not tested * fit: does not converge, or does not estimate correctly for some distributions, this is expected since maximum likelihood does not always work * distributions that have problems for some range of parameters * incorrect random number generator for one discrete distribution Is it useful to put a warning about these issues in stats.distributions in the release? A friend of mine is doing value-at-risk calculations for operational risk analysis, but I am still reluctant to recommend scipy.stats, when it is not clear whether the results are correct or not, although scripting in Python would be much easier than to struggle with some commercial package with a GUI user interface. The problem is that in this application they are using many distributions automatically, and not just a few that are well tested by frequent use. Per Brodtkorb has an improved version of distributions.py with an additional estimation procedure, which we can discuss as an enhancement after 0.7 is out of the way, if there is an interest in this. Josef

4 15

Scipy workflow (and not tools).
by Matthew Brett 16 Mar '09

16 Mar '09

Hi, I've split this off into a new thread because I felt there were two issues in Stefan's original thread. This is in the hope that we can stimulate discussion on the workflow (as opposed to - say - which version control system to use, or which bugtracker). I would be very interested to see if we can come to a consensus on the important discussion of whether to introduce fairly formal code review into the scipy workflow. I've appended the key piece of discussion below. > 2009/2/24 Travis E. Oliphant <oliphant(a)enthought.com>: >> I think the biggest problem has been time and adding too formal of a >> process will just increase the time it takes to get code into SciPy. >> I'm fine with emphasizing documentation and tests as we discuss things >> and we should encourage each other, but I'm not comfortable with >> hard-line statements like the ones being made above. Yes, such things >> are helpful, but they are also expensive and I worry more about what we >> lose in contributions. > > Having so little time means that we cannot be cavalier about adding > broken code to SciPy. Like Matthew mentioned, this becomes an immense > maintenance burden. > >> The quality of what we create should emerge as all interested parties >> critically look at the code that is available in SciPy. > > I agree with that sentiment; and looking critically at code in SciPy > starts with our own patches. > >> Not everyone >> can do that on the same schedule. I'm opposed to trying to force that >> to happen. I very much favor cultivating a culture that wants someone >> to fix the problems in their code. > > Sure, let's be inclusive, but also set a bar. If you make the time to > write a patch, make the time to do it well (it does not take long to > construct a test -- you have to make sure your code works properly > anyhow). > >> But, my favorite workflow is a bit more chaotic, than that. People >> create their own DVCS versions of SciPy using their best judgment and >> publish revisions they consider to be working code. >> >> Branches that are given the thumbs up by 2 people (or 1 on the steering >> committee) get pushed to the main branch. This review happens >> regularly, on IRC channels at regularly scheduled times. > > Two eyes on every piece of code in SciPy, that's all we need. Two > critical eyes that realise the value of tests and documentation. Your > outline above fits in with my view of how this could happen. Best, Matthew

25 96

Changing the container for SciPy's FFTPACK cache?
by Sturla Molden 16 Mar '09

16 Mar '09

SciPy's interface to fftpack is different from NumPy's fftpack_lite. One notable difference is that SciPy's uses a cache in C, whereas NumPy uses a Python dict for the same purpose. SciPy's cache is accessed directly from the C wrappers to FFTPACK. This means that SciPy's FFTs are not threadsafe, and thus the GIL cannot be released. Thus, SciPy has to lock up the interpreter while doing FFTs. Using a Python container with ndarrays as cache would allow the GIL to be released when calling FFTPACK (cf. ticket #1055 for NumPy). I've posted a patch for NumPy yesterday, and I thought of giving SciPy a shot as well (when I can find some spare time). Sturla Molden

1 0

#880: lfilter, variable number of output
by David Cournapeau 15 Mar '09

15 Mar '09

Hi, A user of scipy.signal.lfilter reported a bug in it, which I fixed, but he also raised an interesting point concerning the variable number of output in lfilter. I think always outputting the final values of the delay line is better as well, but we can't change it wo breaking lfilter API. I first thought about adding a lfilter2, but the name is a bit misleading - and I also noticed that filter is not taken (which is the name of this functionality under matlab). What about adding a filter function which is like lfilter, but always return zf ? cheers, David

2 2

Question to Travis: what is rdist about?
by Yaroslav Halchenko 13 Mar '09

13 Mar '09

Travis, I wonder if you get a moment and desire to give a bit of theory/history for the hungry people ;) In a recent thread (a part of it is below the body of this email) dealing with instabilities of rdist Josef asked what is the application domain of rdist distribution... he heard about relation to correlation, I mentioned that it is related to the distribution of a coordinate of points on c-dimensional sphere. But I wonder -- what was the original reason for this distribution to appear? where have you found it, or in other words -- what literature source describes it? thanks to git I found that you introduced it in commit 8ce8603696448c171c186ea2aab158cf34e25441 Author: travo <travo@d6536bca-fef9-0310-8506-e4c0a848fbcf> Date: Fri Nov 22 09:04:46 2002 +0000 Changed statistics module to use clasasses. git-svn-id: http://svn.scipy.org/svn/scipy/trunk@648 d6536bca-fef9-0310-8506-e4c0a848fbcf but I can't figure out if it was really a new distribution or refactored from some other one. Thank you in advance! On Fri, 13 Mar 2009, Yaroslav Halchenko wrote: > > >> Google search for r distribution is pretty useless, and I have not yet > > >> found a reference or an explanation of the rdist and its uses. > > > there was just a single page which I ran to which described rdist and > > > plotted sample pdfs. but can't find it now > > I read somewhere, I don't remember where that rdist is the distribution > > of the correlation coefficient, but without more information that's pretty > > useless > doh! sure it is related... hence the name rdist, since pearsons corr > coeff is abbreviated as 'r' ;) hence rdist ;) > http://en.wikipedia.org/wiki/Correlation_coefficient > says that > The distribution of the correlation coefficient has been examined by R. > A. Fisher[2][3] and A. K. Gayen.[4] > but those are 100 and 50 years old books... not sure if we have them > online to check if they were the one who brought analytic function for > it... > and it seems that it is related to the 'multidimensional' correlation > mentioned in the wikipedia > but it is now clear how sample size "fits into equation"... c seems to > relate to the dimensions of the data... > is it possible to trace back who introduced this lovely piece into > scipy? ;) may be we could ask the author? ;) -- .-. =------------------------------ /v\ ----------------------------= Keep in touch // \\ (yoh(a)|www.)onerussian.com Yaroslav Halchenko /( )\ ICQ#: 60653192 Linux User ^^-^^ [175555]

3 2

Re: [SciPy-dev] [Numpy-discussion] A module for homogeneous transformation matrices, Euler angles and quaternions
by Gareth Elston 11 Mar '09

11 Mar '09

Does anyone know any good internet references for defining and using homogeneous transformation matrices, especially oblique projection matrices? I'm writing some tests for transformations.py and I'm getting unexpected results, quite possibly because I'm making naive assumptions about how to use projection_matrix(). Thanks, Gareth. On Thu, Mar 5, 2009 at 3:28 AM, Jonathan Taylor <jonathan.taylor(a)utoronto.ca> wrote: > Looks cool but a lot of this should be done in an extension module to > make it fast. Perhaps starting this process off as a separate entity > until stability is acheived. I would be tempted to do some of this > using cython. I just wrote found that generating a rotation matrix > from euler angles is about 10x faster when done properly with cython. > > J. > > On Wed, Mar 4, 2009 at 5:10 PM, Gareth Elston > <gareth.elston.floss(a)googlemail.com> wrote: >> I found a nice module for these transforms at >> http://www.lfd.uci.edu/~gohlke/code/transformations.py.html . I've >> been using an older version for some time and thought it might make a >> good addition to numpy/scipy. I made some simple mods to the older >> version to add a couple of functions I needed and to allow it to be >> used with Python 2.4. >> >> The module is pure Python (2.5, with numpy 1.2 imported), includes >> doctests, and is BSD licensed. Here's the first part of the module >> docstring: >> >> """Homogeneous Transformation Matrices and Quaternions. >> >> A library for calculating 4x4 matrices for translating, rotating, mirroring, >> scaling, shearing, projecting, orthogonalizing, and superimposing arrays of >> homogenous coordinates as well as for converting between rotation matrices, >> Euler angles, and quaternions. >> """ >> >> I'd like to see this added to numpy/scipy so I know I've got some >> reading to do (scipy.org/Developer_Zone and the huge scipy-dev >> discussions on Scipy development infrastructure / workflow) to make >> sure it follows the guidelines, but where would people like to see >> this? In numpy? scipy? scikits? elsewhere? >> >> I seem to remember that there was a first draft of a guide for >> developers being written. Are there any links available? >> >> Thanks, >> Gareth.

1 0

build_ext fails to pick up fortran compiler
by Avi Purkayastha 10 Mar '09

10 Mar '09

Hi, When building scipy with a new compiler (pathscale) added to the numpy*distutils*fcompiler list, the build log shows that the build process picks up the fortran compiler for parts of the build where needed except for build_ext. Some parts of the log is listed below for explanation. : customize PathScaleFCompiler customize PathScaleFCompiler using build_clib building 'dfftpack' library compiling Fortran sources Fortran f77 compiler: /opt/pathscale/ice9_native_3.2n_B_sicortex/bin/pathf95 -fixedform -O3 Fortran f90 compiler: /opt/pathscale/ice9_native_3.2n_B_sicortex/bin/pathf95 -O3 Fortran fix compiler: /opt/pathscale/ice9_native_3.2n_B_sicortex/bin/pathf95 -fixedform -O3 creating build/temp.linux-mips64-2.4 creating build/temp.linux-mips64-2.4/scipy creating build/temp.linux-mips64-2.4/scipy/fftpack creating build/temp.linux-mips64-2.4/scipy/fftpack/src creating build/temp.linux-mips64-2.4/scipy/fftpack/src/dfftpack compile options: '-c' pathf95:f77: scipy/fftpack/src/dfftpack/dcosqi.f pathf95:f77: scipy/fftpack/src/dfftpack/dcosqf.f pathf95:f77: scipy/fftpack/src/dfftpack/zfftf.f : and the build continues with success on picking up and building with the pathscale compiler on other pieces that need the fortran compiler until.. : running build_ext customize UnixCCompiler customize UnixCCompiler using build_ext resetting extension 'scipy.integrate._odepack' language from 'c' to 'f77'. resetting extension 'scipy.integrate.vode' language from 'c' to 'f77'. resetting extension 'scipy.lib.blas.fblas' language from 'c' to 'f77'. resetting extension 'scipy.odr.__odrpack' language from 'c' to 'f77'. extending extension 'scipy.sparse.linalg.dsolve._zsuperlu' defined_macros with [('USE_VENDOR_BLAS', 1)] extending extension 'scipy.sparse.linalg.dsolve._dsuperlu' defined_macros with [('USE_VENDOR_BLAS', 1)] extending extension 'scipy.sparse.linalg.dsolve._csuperlu' defined_macros with [('USE_VENDOR_BLAS', 1)] extending extension 'scipy.sparse.linalg.dsolve._ssuperlu' defined_macros with [('USE_VENDOR_BLAS', 1)] customize UnixCCompiler customize UnixCCompiler using build_ext customize PathScaleFCompiler warning: build_ext: f77_compiler=pathscale is not available. : scgcc -shared build/temp.linux- mips64-2.4/scipy/cluster/src/hierarchy_wrap.o build/temp.linux- mips64-2.4/scipy/cluster/src/h ierarchy.o -Lbuild/temp.linux-mips64-2.4 -o build/lib.linux- mips64-2.4/scipy/cluster/_hierarchy_wrap.so building 'scipy.fftpack._fftpack' extension warning: build_ext: extension 'scipy.fftpack._fftpack' has Fortran libraries but no Fortran linker found, using default linker compiling C sources : and ultimately failure because of this reason.. scgcc -shared build/temp.linux- mips64-2.4/scipy/interpolate/src/_fitpackmodule.o -Lbuild/temp.linux- mips64-2.4 -lfitpack -o build/lib.linux-mips64-2.4/scipy/interpolate/_fitpack.so building 'scipy.interpolate.dfitpack' extension error: extension 'scipy.interpolate.dfitpack' has Fortran sources but no Fortran compiler found Any suggestions on why build_ext is failing on picking up the fortran compiler or any work-arounds for this? Thanks Avi

1 0

there is a bug in PyArray_CheckFromAny ?
by Luis Saavedra 10 Mar '09

10 Mar '09

Hi, if o is a "vector" <numpy.ndarray of shape (n,)> and <dtype=int32>, the following lines: po = PyArray_CheckFromAny(o,PyArray_DescrFromType(NPY_INT),0,0, NPY_FORCECAST|NPY_OUT_FARRAY|NPY_ELEMENTSTRIDES,NULL); int *ver = (int *)((PyArrayObject *)po)->data; printf("0: %d\n",ver[0]); printf("1: %d\n",ver[1]); show that: 0: 0 1: -1 regardless of the value (size>2). But if I modify the shape of "o" before the previous line, for example shape = (1,n), works fine! Regards, Luis. PD: In debian squeeze: $ dpkg -s python-numpy |grep Version Version: 1:1.2.1-1 Python-Version: 2.4, 2.5

1 0