Mailman 3 March 2009 - SciPy-Dev

Possible Error in Kendall's Tau (scipy.stats.stats.kendalltau)
by Almer S. Tigelaar March 19, 2009

March 19, 2009

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 … [View More]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 [View Less]

4 34

Status on ubuntu jaunty 64bits.
by Xavier Gnata March 18, 2009

March 18, 2009

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 … [View More] 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 [View Less]

2 2

Warning about remaining issues in stats.distributions ?
by josef.pktd＠gmail.com March 16, 2009

March 16, 2009

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 … [View More]

4 15

Scipy workflow (and not tools).
by Matthew Brett March 16, 2009

March 16, 2009

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 … [View More]

25 96

Changing the container for SciPy's FFTPACK cache?
by Sturla Molden March 16, 2009

March 16, 2009

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 … [View More]

1 0

#880: lfilter, variable number of output
by David Cournapeau March 15, 2009

March 15, 2009

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 … [View More]

2 2

Question to Travis: what is rdist about?
by Yaroslav Halchenko March 13, 2009

March 13, 2009

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? … [View More]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] [View Less]

3 2

Re: [SciPy-dev] [Numpy-discussion] A module for homogeneous transformation matrices, Euler angles and quaternions
by Gareth Elston March 12, 2009

March 12, 2009

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 … [View More]

1 0

build_ext fails to pick up fortran compiler
by Avi Purkayastha March 10, 2009

March 10, 2009

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.… [View More]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 [View Less]

1 0

there is a bug in PyArray_CheckFromAny ?
by Luis Saavedra March 10, 2009

March 10, 2009

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 … [View More]

1 0