Mailman 3 March 2013 - SciPy-User

a routine for fitting a straight line to data
by David Pine March 31, 2013

March 31, 2013

I would like to submit a routine to scipy for performing least squares fitting of a straight line f(x) = ax + b to an (x,y) data set. There are a number of ways of doing this currently using scipy or numpy but all have serious drawbacks. Here is what is currently available, as far as I can tell, and what seem to me to be their drawbacks. 1. numpy.polyfit : a. It is slower than it needs to be. polyfit uses matrix methods that are needed to find best fits to general polynomials (quadratic, cubic, quartic, and higher orders), but matrix methods are overkill when you just want to fit a straight line f(x)= ax + b to data set. A direct approach can yield fits significantly faster. b. polyfit currently does not allow using absolute error estimates for weighting the data; only relative error estimates are currently possible. This can be fixed, but for the moment it's a problem. c. New or inexperienced uses are unlikely to look for a routine to fit straight lines in a routine that is advertised as being for polynomials. This is a more important point than it may seem. Fitting data to a straight line is probably the most common curve fitting task performed, and the only one that many users will ever use. It makes sense to cater to such users by providing them with a routine that does what they want in as clear and straightforward a manner as possible. I am a physics professor and have seen the confusion first hand with a wide spectrum of students who are new to Python. It should not be this hard for them. 2. scipy.linalg.lstsq a. Using linalg.lstsq to fit a straight line is clunky and very slow (on the order of 10 times slower than polyfit, which is already slower than it needs to be). b. While linalg.lstsq can be used to fit data with error estimates (i.e. using weighting), how to do this is far from obvious. It's unlikely that anyone but an expert would figure out how to do it. c. linalg.lstsq requires the use of matrices, which will be unfamiliar to some users. Moreover, it should not be necessary to use matrices when the task at hand only involves one-dimensional arrays. 3. scipy.curve_fit a. This is a nonlinear fitting routine. As such, it searches for the global minimum in the objective function (chi-squared) rather than just calculating where the global minimum is using the analytical expressions for the best fits. It's the wrong method for the problem, although it will work. Questions: What do others in the scientific Python community think about the need for such a routine? Where should routine to fit data to a straight line go? It would seem to me that it should go in the scipy.optimize package, but I wonder what others think. David Pine

5 5

How to use Scipy to solve a Two-Point Boundary Value problem, of a nth system Nonlinear Second-Order Differential Equation
by freeblue March 30, 2013

March 30, 2013

Hello to everyone here, I am trying to find the root of the next equation using Newton method: http://www.wolframalpha.com/input/?i=d%2Fdx%28du%2Fdx%29+%3D+%28-3%2F%28k1%… ^ %281%2F3%29%29+*%28%28du%2Fdx%29^%282%2F3%29%29 which has 2 boundary conditions (u(x=0)=0, u(x=n(max) = m (constant)) ... k1, k2, k3, k4 are arrays, calculated on each x gridpoint Please, do someone knows if it is possible to solve this nth nonlinear second-order differential equation using scipy.integrate.odeint? (nonlinear ODE BVP 1-D. I tried to use Central differencing to "simplify" it but the power number (2/3) makes its hard to use later Newton method. Could be possible then to use scipy? Do you know or see somewhere any relative example ? I am trying to find and write a solution on this the last 3 weeks, so any help will be more than welcome!!!!! thank you in advance, Kas

1 0

special's hidden treasures: rootfinding for distribution cdfs
by josef.pktd＠gmail.com March 29, 2013

March 29, 2013

The fortran code for the distribution functions also has rootfinding build in. We can solve for any value given the other values. For example, solve for the non-centrality parameter for the non-central t-distribution. I never thought I would ever need that. Is this the right function for special.nctdtrinc and friends? I didn't manage to follow the trail. \scipy\special\cdflib\cdftnc.f C SUBROUTINE CDFTNC( WHICH, P, Q, T, DF, PNONC, STATUS, BOUND ) C Cumulative Distribution Function C Non-Central T distribution C C Function C C Calculates any one parameter of the noncentral t distribution give C values for the others. C C Arguments C C WHICH --> Integer indicating which argument C values is to be calculated from the others. C Legal range: 1..3 C iwhich = 1 : Calculate P and Q from T,DF,PNONC C iwhich = 2 : Calculate T from P,Q,DF,PNONC C iwhich = 3 : Calculate DF from P,Q,T C iwhich = 4 : Calculate PNONC from P,Q,DF,T I need iwhich=4 solving for the non-centrality parameter >>> t = stats.nct._ppf(0.1, 9, 10) >>> special.nctdtrinc(9, 0.1, t) 10.000000055648481 confidence interval for non-centrality parameter (AFAIU) >>> special.nctdtrinc(24, [0.975, 0.025], 25) array([ 17.68259, 32.26892]) upper value differs slightly from an R package (which has it's own implementation) cross check >>> nc = special.nctdtrinc(24, 0.025, 25) >>> stats.nct.cdf(25, 24, nc) 0.024999999999859579 a paper refers to a SAS function which seems to do the same '''' following SAS syntax: LowNC_CV = TNONCT(tobs., ν, 1 - α / 2), ''' Background: If I understand this part of some papers correctly, then I can get confidence intervals for effect sizes for statistical tests that are based on the t-distributions with non-central t-distribution under the alternative, with essentially just a call to scipy.special. ------------------ It's nice to find the answer in scipy special when we need to look for something we never knew that it existed. **Scipy special needs your help.** ========================== \scipy\special\add_newdocs.py (maybe my checkout is a bit outdated) <...> add_newdoc("scipy.special", "ncfdtridfd", """ """) add_newdoc("scipy.special", "ncfdtridfn", """ """) add_newdoc("scipy.special", "ncfdtrinc", """ """) add_newdoc("scipy.special", "nctdtr", """ """) add_newdoc("scipy.special", "nctdtridf", """ """) add_newdoc("scipy.special", "nctdtrinc", """ """) add_newdoc("scipy.special", "nctdtrit", """ """) <...> The information is available in the fortran docstrings. A pull request that fills in some blanks is an easy way to become a scipy contributer. Josef

3 3

Find and interpolate zero value
by Florian Lindner March 29, 2013

March 29, 2013

Hello, I have an array like: array([[-1. , 2. ], [-0.5, 4. ], [ 0.5, 6. ]]) first row are coordinates, second row is data. Is there a scipy/numpy way to find zero from the first row and interpolate a data value from the second row? Here I the result would be 5 ( 4 + (6-4)/( 0.5 - (-0.5)) * 0.5 ) How would this best be done with scipy/numpy? Thanks, Florian

2 2

Best way to feed data to Gnuplot?
by Florian Lindner March 29, 2013

March 29, 2013

Hello, I want to utilise gnuplot for my plots. I know about matplotlib, but I don't like it (interface, appearance). What is the best way to feed data from a numpy array to gnuplot? At the moment I try to use it like gnuplot -e "plot '-' ...." and write the data to stdin: 1,2,3 2,6,5 e 10,20,30 Is there a better way? What is the best to print the numpy data in the correct format? There is a Gnuplot.py package, however that seems very outdated, last version is from 2008. According to the docs it still depends on Numeric. Thanks! Florian

5 5

Strange behaviour of scipy.test()
by Nils Wagner March 28, 2013

March 28, 2013

Hi all, I am confused by the current behavior of scipy.test() The number of tests is 0. Any idea ? Nils nwagner@linux-itfv:~> python Python 2.7.2 (default, Aug 19 2011, 20:41:43) [GCC] on linux2 Type "help", "copyright", "credits" or "license" for more information. >>> import scipy >>> scipy.test() Running unit tests for scipy NumPy version 1.8.0.dev-1a816c7 NumPy is installed in /home/nwagner/local/lib64/python2.7/site-packages/numpy SciPy version 0.13.0.dev-9ba2fda SciPy is installed in /home/nwagner/local/lib64/python2.7/site-packages/scipy-0.13.0.dev_9ba2fda-py2.7-linux-x86_64.egg/scipy Python version 2.7.2 (default, Aug 19 2011, 20:41:43) [GCC] nose version 1.1.2 ---------------------------------------------------------------------- Ran 0 tests in 0.055s OK <nose.result.TextTestResult run=0 errors=0 failures=0>

2 1

Sequence not work
by El suisse March 27, 2013

March 27, 2013

Hi I have to represent the following sequence: [image: z[2]=2] [image: z[n+1]=2^{n-1/2}\sqrt{1-\sqrt{1-4^{1-n}z[n]^{2}}] [image: n=2,3,4....] and my code is as follows: #!/usr/bin/env python import matplotlib.pyplot as plt import numpy as np num_muestras = 100 z = np.zeros(num_muestras) z[2] = 2.0 #En Python se cuenta desde el cero for n in range(2,num_muestras-1): z[n+1] = np.power(2,n-(1/2)) * (np.sqrt(1- np.sqrt(1-np.power(4,1-n) * np.power(z[n],2))) ) fig1 = plt.figure() ax = plt.gca() plt.ylabel('$z[n+1]=2^{n-1/2}\sqrt{1-\sqrt{1-4^{1-n}z[n]^{2}}}$') plt.xlabel('$n$') plt.grid() ax.plot(z) plt.show() but does not work, z should tend to pi. might be wrong? thanks in advance

3 2

Size of sparse boolean coo_matrix unexpectedly tiny
by DoanVietTrungAtGmail March 26, 2013

March 26, 2013

I created an NxN (N = 10^6) sparse boolean coo_matrix A and randomly filled it with 200,000 Trues. sys.getsizeof(A) said A took 32 bytes - a thousand times smaller than I expected. I am new to Scipy, please tell me what's going on, thanks. >>> rows = [random.randint(0,999999) for j in xrange(200000)] >>> cols = [random.randint(0,999999) for j in xrange(200000)] >>> vals = ones_like(rows) >>> A = sparse.coo_matrix((vals,(rows,cols)), shape = (int(1E6), int(1E6)), dtype = bool) >>> len(vals) 200000 >>> print sys.getsizeof(A) 32 >>>

1 0

constrained optimization: question about Jacobian
by Andreas Hilboll March 26, 2013

March 26, 2013

I want to perform constrained optimization of the function def F(X, data): # x1, x2, x3, t2 are scalars, data is a constant np.array x1, x2, x3, t2 = X T = np.arange(X.size) t1, t3 = T[0], T[-1] Xbefore = x1 + (T - t1) * (x2 - x1) / (t2 - t1) Xafter = x2 + (T - t2) * (x3 - x2) / (t3 - t2) Xbreak = np.where(T <= t2, Xbefore, Xafter) return ((Xbreak - data)**2).sum() where the parameter value t2 must be within 0 <= t2 < data.size. According to the tutorial (http://docs.scipy.org/doc/scipy/reference/tutorial/optimize.html#constraine…), I need to define the Jacobian. However, I'm unsure of how to define the derivative with regards to t2. Any ideas are greatly appreciated :) Cheers, A.

2 2

Re: [SciPy-User] cobyla
by The Helmbolds March 26, 2013

March 26, 2013

On my system (Windows 7, Python 2.7.x and IDLE, latest SciPy), I observe the following behavior with fmin_cobyla and minimize's COBYLA method. Case 1: When run either in the IDLE interactive shell or within an enclosing Python program: 1.1. The fmin_cobyla function never returns the Results dictionary, and never displays it to Python's stdout. This is true regardless of the function call's disp setting. 1.2. The 'minimize' function always returns the Results dictionary but never displays it to Python's stdout. Again, this is true regardless of the function call's disp setting. Case 2: When run interactively in Window's Command Prompt box: 2.1 The fmin_cobyla function never returns the Result dictionary, regardless of the function call's disp setting. Setting disp to True or False either displays the Results dictionary in the command box or not (respectively). I don't think the Results dictionary gets to the command box via stdout. 2.2 The 'minimize' function always returns the Result dictionary, regardless of the function call's disp setting. Setting disp to True or False either displays the Results dictionary in the command box or not (respectively). I don't think the Results dictionary gets to the command box via stdout. My thanks to all who helped clarify this situation. Bob H

3 4