Mailman 3 October 2012 - SciPy-User

stats.ranksums vs. stats.mannwhitneyu
by Nils Kölling Oct. 11, 2012

Oct. 11, 2012

I am trying to perform a Mann-Whitney U (AKA rank sum) test using Scipy. My data consists of around 30 samples in total with ties, so I get anything between 1:29 .. 15:15 .. 29:1 samples per group. As far as I can see there are two options: scipy.stats.ranksums: Does not handle ties, equivalent to R's wilcox.test with exact=False and correct=False scipy.stats.mannwhitneyu: Handles ties, equivalent to R's wilcox.test with exact=False and correct=use_continuity So at first glance the MWU function would seem to be the better choice, except the docs explicitly state that it should not be used with less than 20 samples per group. So what is the best function to use in this case? What kind of biases will I get when I use the mannwhitneyu function with less than 20 samples? And what sort of problems do ties cause with ranksums? Cheers Nils

3 5

special: beta or gamma
by josef.pktd＠gmail.com Oct. 10, 2012

Oct. 10, 2012

scipy.stats has a bit of a mixed usage in most cases we use gammaln to get higher precision, but maybe betaln would be better >>> factorial(30) / factorial(5) / factorial(25) 142506.0 >>> 1./special.beta(6, 26) / 31 142505.99999999997 >>> special.gamma(31) / special.gamma(6) / special.gamma(26) 142506.0 >>> np.log(factorial(30) / factorial(5) / factorial(25)) 11.867139383067599 >>> -special.betaln(6, 26) - np.log( 31) 11.867139383067599 >>> special.gammaln(31) - special.gammaln(6) - special.gammaln(26) 11.867139383067581 or maybe there is no difference >>> n1, n2 = 5000, 4000 >>> -special.betaln(n1+1, n2+1) - np.log(n1 + n2 + 1) 6177.8820911143594 >>> special.gammaln(n1+n2+1) - special.gammaln(n1+1) - special.gammaln(n2+1) 6177.8820911143594 >>> n1, n2 = 50000, 4000 >>> special.gammaln(n1+n2+1) - special.gammaln(n1+1) - special.gammaln(n2+1) 14253.783294794481 >>> -special.betaln(n1+1, n2+1) - np.log(n1 + n2 + 1) 14253.78329479451 some ``special`` notes while finding out how many permutations we have. Josef

1 0

CFP: Scientific Software Days, Dec 17, Austin, TX
by Scientific Software Days Oct. 10, 2012

Oct. 10, 2012

CALL FOR PARTICIPATION: 6th ANNUAL SCIENTIFIC SOFTWARE DAYS Austin TX December 17 2012 Conference Details and Talk Submission at http://scisoftdays.org/meetings/2012/ Please email questions to softwareday(a)tacc.utexas.edu Hosted by the Texas Advanced Computing Center and the Jackson School of Geosciences, University of Texas at Austin. Scientists use software for their research. Some of them also develop computational software as part of their research. Scientific Software Day is an ongoing meeting of users and producers of scientific software, with presentations by scientific software tool makers and the users of their tools. The objective is to build cross-disciplinary community and skills in the diverse set of users and developers of scientific software, both academic and industrial. Most groups that use supercomputing cope with their scientific software environment in isolation, not always relying on prepackaged “canned” solutions. Many successful lines of research and development are achieved, but many times less than optimal paths are taken, simply because computing is done by people stretched between computational skills and skills in the relevant science and engineering specialties. Available tools and methods are not always known to the people who need them, and time pressure makes it hard to make the best use of the tools available. Support staff at supercomputing centers is stretched and is best at responding to specific issues rather than offering broad support. We seek to build a community to address these needs. The Scientific Software Day at UT Austin is intended to nucleate that community. If you are involved in any end use or development of scientific software, you can benefit from and contribute to this goal. This is, therefore, a somewhat unusual call for presentations. Ideal presentations for Scientific Software Days are of two types: 1) presentations of generic tools that can be used in scientific software development and deployment 2) presentations of specific work, focusing on experience in developing scientific software, workflows, and tool chains. We are especially seeking presentations of the second type. We would appreciate a brief introduction to your work intended for a general scientific audience, and then a focus on your workflow or any particular aspect of it that presented particular challenges or required original solutions. The target audience will be a broad selection of the scientific and engineering communities with a particular interest in supercomputing. Let’s get to know each other and learn from one another. Andy R. Terrel, Ph.D. Scientific Software Days Organizer Texas Advanced Computing Center University of Texas at Austin aterrel(a)tacc.utexas.edu

1 0

linalg.iterative problem!
by SungHwan Choi Oct. 10, 2012

Oct. 10, 2012

Hi I have a trouble with iterative linear equation solver in sparse linalg >>> A=np.random.rand(100,100) >>> b=np.random.rand(100,1) >>> x1,info=bicg(A,b) >>> print info 0 >>> x,info=bicg(A,b,x0=x1) >>> print info 0 Wheather setting guessing value or not, solutions should be the same if both calculations were converged. But x is ridiculous values >>> x1 array([ -1.23619915e+00, 2.97586147e+00, -6.14648472e-01, -1.60393960e+00, 1.27603461e+00, -6.71559072e-01, -4.57949314e-04, -2.77060765e-01, -4.17339328e-01, 4.58723928e-01, 1.25520390e+00, -6.43505834e-01, -2.50714182e+00, 1.64801812e+00, 3.22769225e-01, -3.47764545e+00, -2.64155124e+00, 8.02161189e-01, -2.04397410e-02, 1.78176386e+00, -2.30534938e+00, 2.03747784e-01, 4.25741370e-01, 8.92017739e-02, 7.92235549e-01, 2.05296800e+00, -5.07849138e-01, 2.39548767e+00, -6.75288110e-01, 5.40248358e-01, -1.22652305e+00, 1.24128988e+00, 3.55832137e-01, -4.94905114e-01, 1.89255642e+00, 1.13032169e+00, -1.13126641e+00, -1.24107851e+00, -3.50610928e-01, 1.51242380e+00, 5.93313109e-02, -1.65542281e+00, -1.31457525e+00, -2.07950912e+00, 2.03842426e+00, -1.25129931e+00, -1.18204676e+00, -2.84095828e-01, 1.50420723e+00, -1.86947284e+00, 3.82634122e-01, 1.59583715e+00, 2.38088734e+00, -1.94456801e+00, -3.91679300e+00, -5.82275859e-01, 6.37373111e-01, 1.50117747e+00, 3.16166509e-01, -4.80709301e-01, 2.44748482e-01, 8.46311114e-01, -3.50561001e-01, 1.17040825e+00, 8.48462084e-01, 2.26995940e+00, -4.02400162e-01, 8.15964837e-02, -1.17082091e-01, 5.76520318e-01, 2.68571769e+00, -8.24618021e-01, 1.70237224e+00, -9.51878209e-01, -1.79056788e+00, 1.28023233e+00, -3.06323112e+00, 1.36928031e+00, -5.32667426e-01, -8.76808999e-01, 5.04791986e+00, 1.02573111e+00, 2.91480759e-01, -1.65205205e+00, 2.01570733e+00, -8.58303160e-01, 1.00844953e+00, -1.42026281e+00, -1.35743978e+00, -6.98618293e-01, -9.64603408e-01, 4.94354222e-01, -1.56639931e+00, -1.00424343e+00, 1.34539380e+00, 8.34746938e-01, -1.42944790e-02, 4.11728888e-02, 8.48928870e-01, -3.81714583e-01]) >>> x array([ 5.54091698e-06, 1.01097308e-06, -1.90554231e-07, 4.72137778e-07, -1.59217304e-06, -2.63521296e-06, 1.82539010e-06, -3.29784996e-06, 2.11241995e-06, 3.34528259e-06, -1.48936133e-06, 4.64833156e-06, 2.72517397e-06, -8.68280493e-07, -1.48461475e-06, 1.31078987e-06, 1.96827837e-06, 2.43522800e-06, -1.81519616e-09, 7.36595257e-07, -1.68678301e-06, -2.36489475e-06, -9.48767026e-08, -4.19287423e-07, -1.94382913e-06, -2.85541661e-06, -2.22431928e-06, 5.69426787e-07, -3.20549054e-06, -4.28991209e-06, -2.66204912e-06, -5.41291369e-07, 1.40179165e-07, -7.73036341e-08, -2.62207353e-06, 3.04217252e-07, 7.58099103e-06, 1.64647208e-07, -2.07367685e-06, 1.07293388e-06, -2.64252934e-06, 8.43832882e-07, -2.09558797e-06, -2.38424059e-06, -2.01101471e-06, 1.14992748e-06, 1.75975671e-06, 3.47029359e-06, -1.73474476e-06, 1.63282775e-06, 2.14847352e-06, -1.06630511e-07, -3.71185399e-07, -6.19298483e-07, -4.22283992e-07, 2.87057463e-06, 2.50493018e-07, 2.38959629e-07, 1.09429464e-06, 2.78931839e-06, 1.04950522e-06, 1.92574749e-06, 2.16166697e-06, -1.49381992e-07, 2.57534472e-06, -1.80238481e-06, -4.48006258e-07, -4.82004956e-06, -2.90858804e-06, 2.36872252e-06, -5.82462798e-06, 2.28721650e-06, 4.98778955e-06, 7.42277728e-07, 4.79308235e-06, 3.32154978e-06, 2.01826593e-06, 1.70133451e-06, 1.04876888e-06, -1.66519455e-06, -2.29374493e-06, 2.85916887e-06, -2.92097942e-06, 4.34734275e-07, 2.09635331e-06, -1.21218109e-06, -2.15483189e-06, -2.62789759e-06, 5.97557810e-06, -8.16033223e-07, 5.59003423e-07, -1.27845573e-06, -2.81987257e-06, -2.99137673e-06, -1.68537057e-06, -1.23698610e-06, -1.26839543e-06, -6.37207587e-08, 3.99191552e-07, -1.97820929e-06]) When I set x0, I always get some very small value as solution but I don't know why it give us original solution Please, help me if you have a piece of knowhow to this phenomena Sincerely Sunghwan

1 0

Just curious
by The Helmbolds Oct. 9, 2012

Oct. 9, 2012

Just curious as to why 'anneal' is one of the allowable 'minimize' method options, while 'brute' is not. Bob H

1 0

numpy.piecewise doesn't work with lists, only ndarrays
by Juan Luis Cano Rodríguez Oct. 9, 2012

Oct. 9, 2012

I have noticed this behaviour of numpy.piecewise: In [1]: import numpy as np In [2]: q = [1, 2, 3, 4, 5, 6] In [3]: np.piecewise(q, [q < 3, 3 <= q], [-1, 1]) Out[3]: array([ 1, -1, 0, 0, 0, 0]) In [4]: q = np.array(q) In [5]: np.piecewise(q, [q < 3, 3 <= q], [-1, 1]) Out[5]: array([-1, -1, 1, 1, 1, 1]) Maybe the function should work the same both with lists and arrays? Should I file a bug?

4 5

Thoughts on SciPy Conference
by Andy Ray Terrel Oct. 9, 2012

Oct. 9, 2012

Hello all, I wrote up some thoughts on the SciPy conference. I would appreciate any feedback as I'm chatting with some Enthought employees about next year. http://andy.terrel.us/blog/2012/10/08/thoughts-on-the-scipy-conference/ -- Andy

1 0

Re: [SciPy-User] new scipy website
by Thomas Kluyver Oct. 8, 2012

Oct. 8, 2012

Hi Joris, Thanks for all the thoughts. I've CCed the scipy-user list in this response, so other people can pitch in as well. On 8 October 2012 22:49, Joris Van den Bossche <jorisvandenbossche(a)gmail.com> wrote: > Dear Thomas, > > As you said on the numfocus list, there are a lot out there with interest in > getting involved in the scipy website, and maybe I am someone like that. > > I am a PhD Student at Ghent University, and an enthousiast scientific python > user, also following the lists, but never really getting involved (apart > from sporadic bug reporting). And I also don't have that much of a time, but > maybe every little bit helps. > > But at the moment, that is my feeling at least with the new scipy website, > it is not very clear what has to be done, some kind of 'roadmap' for the > website, where interested people can get involved. > There were some discussions on the mailing lists, like the one in end of > July (http://thread.gmane.org/gmane.comp.python.scientific.devel/16796) , > but never really concrete. There was also no activity the last 5 months at > the scipy.org-new repository. I don't want to leap in and present a whole roadmap straight away - as we saw with the packages question, people have a wide range of opinions on how best to achieve the same basic aim, and I'm sure opinions on the website will be similarly diverse. I've largely used up my disagreement quota for a couple of weeks. ;-) I'm assuming that we want to work roughly within the framework of the existing 'new' site, which is a set of static pages built by Sphinx (as is ipython.org). Although many of us are familiar with Sphinx, there is a case for moving to a tool more suited to websites than documentation. If someone wanted to do some research and present the case for such a change, I think we'd seriously consider it. > At the numfocus list, you mentioned some first to do points: > >> - Update the front page to promote the Scipy stack we've agreed on: >> adding Pandas & Sympy, rearranging the current distinction of 'core >> projects' vs. 'related projects'. >> - A page describing the Scipy stack specification. >> - A new separate page (or pages) about scipy-the-package. >> - Some general design work wouldn't go amiss - do we need the >> breadcrumb bar? Can we improve top level navigation? One more thing people can get involved with straight away. On the new installation page I've done, there's a section for installing from Linux distro packages. I've filled it in for Ubuntu & Debian. Could users of other major distributions provide a short entry for each, with the command to install all the packages, the first distro version which meets the specification, and any other relevant info? Here's the section for Ubuntu & Debian: https://github.com/scipy/scipy.org-new/pull/3/files#L3R29 > Had you already an idea on how you wanted to proceed? > I was thinking that maybe some kind of discussion document, for discussing > an outline of the design and structure of the site we want in the future > could be helpfull. To engage others in the discussion. For later, to see > what still has to be done. Also, to really think about the structure (which > pages we want, ...) we envision in the future. Not necessarily directly > achievable from the start, but to have a goal to work to, and to avoid that > each time somebody wants to add something adds a page and you end op with a > mess like the scipy.org today (at least, I find it a not very clear and > structured site). I was thinking of a google docs, but maybe something > similar is possible on github. I think such dynamic document is better > suited for this that a discussion on a mailing list. That does sound like a good idea. Github has wiki pages, and of course the current scipy.org site is a wiki. I'm also happy to use a Google Doc, or one of the successors to Google Wave, like Rizzoma. Since you had the idea, can you pick a platform, and sketch out some sort of overview, so that we can start discussing and filling in details? > Apart from the points you raised above, I was also thinking about the > following: > - What do we want with a lot of material on the scipy website right now, eg > cookbook, topical software? Bring it over (but then they need an update I > think)? Do we want something like modern version of the cookbook > (cross-project examples, maybe integrated with notebooks), or should it be > replaced by Scipy central? For now, much of the content will probably stay on the existing server, perhaps at a subdomain, because updating it would take person-hours from more important jobs. Longer term, I hope we can get it moved to new homes as appropriate. Scipy-central unfortunately also seems somewhat dead. We got in touch with the developer recently, and he essentially offered the Django codebase to anyone who has the time to work on it. There's also work going on with nbviewer.ipython.org, which could be the base of a new example-code sharing site. > - Include some specific documentation/tutorial at the scipy stack level? > (maybe based on https://github.com/scipy-lectures/scipy-lecture-notes) Yes, longer term, I hope that we'll develop much more connected documentation around the scipy stack, including tutorials and howtos. Best wishes, Thomas > Sorry for the long mail. It were just some ideas, see what you can do with > it. > > Regards, > Joris

1 0

Scipy stack: standard packages (poll)
by Kyle Mandli Oct. 8, 2012

Oct. 8, 2012

We were having a discussion on this topic on Scipy-Dev where Pauli Virtanen mentioned that there had been effort to implement a plan similar to what iPython and matplotlib does with their documentation and website. The conversation has dropped off from there but I think there is enough there already to get started. The best information on our status now was from Pauli a month ago: http://mail.scipy.org/pipermail/scipy-dev/2012-August/017916.html Kyle > Hey Thomas, > >On Mon, Oct 8, 2012 at 7:21 AM, Thomas Kluyver <takowl(a)gmail.com> wrote: >> >> Thanks, Fernando. I think the next step is to start reworking the >> 'new' scipy.org website (http://scipy.github.com/ ) to focus on the >> stack, rather than scipy-the-package. I've just kicked that off with a >> pull request replacing the 'download' page with an 'install' page: > >Kyle Mandli had already been pushing hard on this front, it would be >great if the two efforts could play off each other, as he'd spent a >good amount of time already on this idea... > >f

1 0

SciPy ODEINT Problem
by Lofgren, Eric Oct. 8, 2012

Oct. 8, 2012

I've been working on a set of ordinary differential equations for an epidemic model. Normally odeint works swimmingly for this kind of thing, though I'll admit this one is a touch more complex than most of the ones I've thrown at it. I'm getting the following error message when I try to run the code below: lsoda-- at t (=r1) and step size h (=r2), the corrector convergence failed repeatedly or with abs(h) = hmin In above, R1 = 0.7016749763132E+04 R2 = 0.1954514552051E-06 Repeated convergence failures (perhaps bad Jacobian or tolerances). Having used the atol and rtol options to relax the tolerances by quite a bit, the error then converts to: RuntimeWarning: overflow encountered in double_scalars This feels to me like a coding error at that point, rather than some issue with the solver itself, but I've not managed to find anything in particular. Any ideas? Thanks, Eric # Imports import numpy as np from pylab import * import scipy.integrate as spi # Initial Population States - Model is in Individuals Us0 = 2.0 H0 = 1.0 Up0 = 4.0 Cp0 = 4.0 Ua0 = 4.0 Ca0 = 4.0 D0 = 0.0 PopIn = (Us0, H0, Up0, Cp0, Ua0, Ca0, D0) # Model parameters # Time is currently in MINUTES N = Us0 + H0 + Up0 + Cp0 + Ua0 + Ca0 + D0 M = Up0 + Cp0 + Ua0 + Ca0 + D0 n_contacts = (3.0)*(1/20.0) p_contacts = nurse_contacts/N rho_p = p_contacts rho_d = p_contacts rho_a = p_contacts sigma_p = 0.05 sigma_d = 0.25 sigma_a = 0.05 omega = 14400 alpha = 0.25 psi_p = 0.10 psi_a = 0.10 mu_p = 0.0 mu_a = 0.0 theta_p = 1.0/10080.0 theta_a = 1.0/10080.0 nu_cp = 0.07 nu_ca = 0.07 nu_d = 0.0 nu_up = 0.43 nu_ua = 0.43 kappa = 0.10 tau = 3.50 iota = 1/20.0 * 0.60 * 1.00 zeta = 0.2785 * (1.0/17280.0) # Pr(death) and time until death gamma = (1.0-0.2785) * (1.0/12902.4) # Pr(discharge) and time until discharge theta = theta_p + theta_a # ODE Fit and Graphing Parameters t_end = 144000 t_start = 1 t_step = 0.1 t_interval = np.arange(t_start, t_end, t_step) # The actual model running part def eq_system(PopIn,t): #Creating an array of equations Eqs= np.zeros((7)) Eqs[0] = ((iota*PopIn[1]) - (rho_p*sigma_p*PopIn[3]*(PopIn[0]/N)) - (rho_d*sigma_d*PopIn[6]*(PopIn[0]/N)) - (rho_a*sigma_a*PopIn[5]*(PopIn[0]/N))) Eqs[1] = ((rho_p*sigma_p*PopIn[3]*(PopIn[0]/N)) + (rho_d*sigma_d*PopIn[6]*(PopIn[0]/N)) + (rho_a*sigma_a*PopIn[5]*(PopIn[0]/N)) - (iota*PopIn[1])) Eqs[2] = (((1/omega)*PopIn[4]) - alpha*PopIn[2] - (rho_p*psi_p*PopIn[2]*(PopIn[1]/N)) - (mu_p*sigma_p*PopIn[2]*(PopIn[3]/N)) - (mu_a*sigma_a*PopIn[2]*(PopIn[5]/N)) - theta_p*PopIn[2] + nu_up*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6]))) Eqs[3] = (alpha*PopIn[2] - ((1/omega)*PopIn[4]) - (rho_a*psi_a*PopIn[4]*(PopIn[1]/N)) - (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) - (mu_a*sigma_a*PopIn[4]*(PopIn[5]/N)) - theta_a*PopIn[4] + nu_ua*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6]))) Eqs[4] = (((1/omega)*PopIn[5]) + (rho_p*psi_p*PopIn[2]*(PopIn[1]/N)) + (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) + (mu_a*sigma_a*PopIn[2]*(PopIn[5]/N)) - alpha*PopIn[3] - kappa*PopIn[3] - theta_p*PopIn[3] + nu_cp*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6]))) Eqs[5] = (alpha*PopIn[3] + (rho_a*psi_a*PopIn[4]*(PopIn[1]/N)) + (mu_p*sigma_p*PopIn[4]*(PopIn[3]/N)) + (mu_a*sigma_a*PopIn[4]*(PopIn[5]/N)) - ((1/omega)*PopIn[5]) - kappa*tau*PopIn[5] - theta_a*PopIn[5] + nu_ca*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6]))) Eqs[6] = (kappa*PopIn[3] + kappa*tau*PopIn[5] - gamma*PopIn[6] - zeta*PopIn[6] + nu_d*((theta*M)+(zeta*PopIn[6])+(gamma*PopIn[6]))) return Eqs # Model Solver model = spi.odeint(eq_system, PopIn, t_interval)

2 1