Mailman 3 March 2010 - SciPy-User

calculating intersection points in probability curves
by Tim Michelsen 03 Mar '10

03 Mar '10

Hello, I am performing a curve fit and also showing a distribution histogram envelope in my plot. In order to help the reader to evaluate the result I would like to draw certain boundaries (vertical and horizontal lines). While I am aware on how to draw such lines with matplotlib, I would like to know whether there are some functions in scipy or matplotlib which help me to retrieve the coordinates a) at which two curves intersect b) at which a distribution reaches a certain value? Example: How do I get the y-axis value which is reached by the green curve in http://matplotlib.sourceforge.net/_images/histogram_demo_extended_021.png at a x-axis value of in 175? I could proably use a solver from numpy like http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.solve.html… but if I plot a distribution, the equation of the envelove is unknown at the first place. I'd appreciate your help or pointers to examples. Thanks a lot in advance, Timmie

2 2

scipy.interpolate.rbf sensitive to input noise ?
by denis 03 Mar '10

03 Mar '10

A followup to the 15feb thread "creating a 3D surface plot from collected data": the doc for scipy.interpolate.rbf says that 3 of the 7 radial functions use a scale paramater `epsilon`, 4 do not. Odd -- how can some be scale-free ? Running rbf on 100 1d random.uniform input points (after changing the linalg.solve in rbf.py to lstsq) shows that all but linear and gaussian are noisy, giving interpolants way outside the input range: N: 100 ngrid: 50 sigma: 0.1 # min max max |delta| Rbf -1.0 1.0 0.5 gaussian -1.0 1.2 0.3 linear -1.4 2.1 2.5 thin-plate -1.0 2.2 3.0 inverse multiquadric -1.0 2.4 3.2 multiquadric -2.1 12.8 7.4 quintic -10.6 7.0 11.7 cubic Should Rbf be moved off to noisymethods/... until experts can look it over, or have I done something stupid ? (A back-of-the-envelope method for choosing epsilon aka r0 would be nice, and a rough noise amplification model would be nice too.) Also, the doc should have a link to http://www.scipy.org/Cookbook/Matplotlib/Gridding_irregularly_spaced_data and should say that the initial rbf() is N^2. cheers -- denis <pre><code> # from http://scipy.org/Cookbook/RadialBasisFunctions import sys import numpy as np from rbf import Rbf # lstsq # from scipy.interpolate import Rbf Rbfuncs = ( 'gaussian', 'linear', 'thin-plate', 'inverse multiquadric', 'multiquadric', 'quintic', 'cubic', ) N = 100 ngrid = 50 sigma = .1 exec "\n".join( sys.argv[1:] ) # N= np.random.seed(1) np.set_printoptions( 2, threshold=50, edgeitems=5, suppress=True ) # .2f x = np.sort( np.random.uniform( 0, 10, N+1 )) y = np.sin(x) + np.random.normal( 0, sigma, x.shape ) xi = np.linspace( 0, 10, ngrid+1 ) print "N: %d ngrid: %d sigma: %.2g" % (N, ngrid, sigma) print "# min max max |delta| Rbf" for func in Rbfuncs: rbf = Rbf( x, y, function=func ) fi = rbf(xi) maxdiff = np.amax( np.abs( np.diff( fi ))) print "%5.1f %4.1f %4.1f %s" % ( fi.min(), fi.max(), maxdiff, func ) </pre></code>

3 16

Scipy's probplot compared to R's qqplot
by PHobson＠Geosyntec.com 03 Mar '10

03 Mar '10

Hey folks, I've taken more of an interest in statistics and Scipy lately and decided to compare the scipy.stats.probplot() function to R's qqplot(). For a given dataset, the results are slightly different. Here's a link to the script I wrote to do the comparison. http://dpaste.com/167464/ Basically, it does the following: -Uses numpy to generate some fake, noramlly distributed data -Uses both R and Scipy to compute the values needed for quantile/probability plot -Computes linear regressions on the quantile data with both R and Scipy. -prints some output to compare the two My initial conclusions: 1) R's lm(y~x) and scipy.stats.linregress(x,y) yield the same slope and intercept of a linear model. (good) 2) R and Scipy compute the quantiles of a dataset in slightly different manners (??) Any clue as to why the discrepancy in #2 occurs? Would you consider it a big deal? I'm using: Python v2.6.2 (XP) and v2.6.4 (Karmic and Snow Leopard) Scipy v0.7.1 Numpy v1.4.0 R v2.10.0 Rpy2 v2.0.8 Thanks, -Paul H.

3 5

segfaults with weave.inline
by pitam 03 Mar '10

03 Mar '10

I did some more research on this. My code uses blitz array, where the original array was a numeric array of numpy. This is where the fault is at. I need to convert the blitz array back to a numeric numpy array. How do I do this, any Ideas?? Original post: http://mail.scipy.org/pipermail/scipy-user/2010-March/024478.html Follow-ups: http://mail.scipy.org/pipermail/scipy-user/2010-March/024498.html Pitam Mitra Deepspace Group UCSB Physics

2 1

Re: [SciPy-User] helmert implementation (7 parameters - geometric transformation)
by Massimo Di Stefano 03 Mar '10

03 Mar '10

Hi, from here : http://docs.scipy.org/doc/numpy/reference/generated/numpy.linalg.cholesky.h… seems it is applicable only on a square matrix. is that true ? reading this pdf : http://geomatica.como.polimi.it/corsi/misure_geodetiche/MG_03SistRif2.pdf scroll to the page : 15 i can see it isn't a common system : Y = AX but it is : Y = AX + B + V the points are : p1: 4402553.569 727053.737 4542823.332 4399375.518 703845.639 4549214.880 4412911.336 701094.214 4536517.955 p2 : 4402553.334 727053.937 4542823.474 4399375.347 703845.876 4549215.105 4412911.150 701094.435 4536518.139 tring to rewrite the code i did : L = np.loadtxt(p1) A = np.zeros((3*L.shape[0],7),float) A[ ::3, 0] = 1.0 A[1::3, 1] = 1.0 A[2::3, 2] = 1.0 A[1::3, 3] = L[:,2] A[2::3, 3] = -L[:,1] A[ ::3, 4] = -L[:,2] A[2::3, 4] = L[:,0] A[ ::3, 5] = L[:,1] A[1::3, 5] = -L[:,0] A[ ::3, 6] = L[:,0] A[1::3, 6] = L[:,1] A[2::3, 6] = L[:,2] G = np.loadtxt(p2) V = np.zeros((3*G.shape[0],1),float) V[ ::3, 0] = G[:,0] - L[:,0] V[1::3, 0] = G[:,1] - L[:,1] V[2::3, 0] = G[:,2] - L[:,2] i used the code in the previouse mail : N = np.dot(A.T.conj(), A) T = np.dot(A.T.conj(), Y) C = np.dot(linalg.inv(N), T) to solve a system like : Y = AX how can i change the code to solve : Y = AX + B + V instead ? the results in the dpf are : x0 -9.256 m y0 -23.701 m z0 16.792 m Rx -0.0001990982 rad Ry 0.0001778762 rad Rz 0.00015 rad μ 0.00000046 any suggestion, also on how to apply a decomposition can give me a great help! Il giorno 02/mar/2010, alle ore 19.00, scipy-user-request(a)scipy.org ha scritto: > Hi, > > I see your coordinates are very big numbers. These can lead to an > ill-conditioned (AtA) matrix which can lead to a wrong inverse. Try > using cholesky decomposition instead of using linalg.inv. The method > is outlined here http://en.wikipedia.org/wiki/Cholesky_decomposition > and numpy has routines for doing just this. The other option is to > subtract the centre-of-mass from your coordinates i.e. P1 -> P1x - > P1xave, P1y - P1yave etc. Do the same for P2. Calculate the > transformation again. The shifts should now be zero but your rotations > and scale should be statistically independent and correct. The > difference between the two coordinate systems' centres of mass gives > the shifts. Transformations far away from the origin rarely behaves > well. > > Good luck.

3 2

scipy.io.wavfile.read
by Nils Wagner 03 Mar '10

03 Mar '10

Hi all, I tried to read a *.wav file using scipy.io.wavfile.read from scipy.io import wavfile from pylab import plot, show samplerate, data = wavfile.read('BagPipes.wav') It failed with /home/nwagner/local/lib64/python2.6/site-packages/pytz/__init__.py:29: UserWarning: Module pkg_resources was already imported from /home/nwagner/local/lib64/python2.6/site-packages/setuptools-0.6c9-py2.6.egg/pkg_resources.py, but /home/nwagner/local/lib64/python2.6/site-packages/distribute-0.6.4-py2.6.egg is being added to sys.path from pkg_resources import resource_stream /home/nwagner/local/lib64/python2.6/site-packages/pytz/__init__.py:29: UserWarning: Module site was already imported from /usr/lib64/python2.6/site.pyc, but /home/nwagner/local/lib64/python2.6/site-packages/distribute-0.6.4-py2.6.egg is being added to sys.path from pkg_resources import resource_stream /home/nwagner/local/lib64/python2.6/site-packages/scipy/io/wavfile.py:20: WavFileWarning: Unfamiliar format bytes warnings.warn("Unfamiliar format bytes", WavFileWarning) /home/nwagner/local/lib64/python2.6/site-packages/scipy/io/wavfile.py:92: WavFileWarning: chunk not understood warnings.warn("chunk not understood", WavFileWarning) Traceback (most recent call last): File "read_wav.py", line 9, in <module> samplerate, data = wavfile.read('BagPipes.wav') File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/io/wavfile.py", line 90, in read data = _read_data_chunk(fid, noc, bits) File "/home/nwagner/local/lib64/python2.6/site-packages/scipy/io/wavfile.py", line 43, in _read_data_chunk data = numpy.fromfile(fid, dtype=dtype, count=size//bytes) ZeroDivisionError: integer division or modulo by zero Any idea ? The example file is available at http://www.ilovewavs.com/Effects/Music/Music.htm Cheers, Nils

2 2

segfaults with weave.inline
by pitam 02 Mar '10

02 Mar '10

Hi, I did some more testing on this. The zone_vector() had a wrong comment on it. Its really small, it is 8*(number of elements in the star shaped pattern - that is zone_outer_radius - zone_inner_radius) I added some extra lines to see where the thing breaks - it does the median filter to the full image. After the last pixel, it breaks. This is about the same time where the C code should end. Hence, my guess is there is a memory leak. IS there any way to check this? Pitam Mitra Deepspace Group, UCSB Physics.

1 0

splmake gone from documentation
by Jason Grout 02 Mar '10

02 Mar '10

I'm teaching a numerical analysis class and we are talking about different kinds of splines. I notice in some old scipy documentation that splmake is a nice function for making different kinds of splines (i.e., clamped, free, not_a_knot, etc.). However, I can't find splmake in the reference documentation at http://docs.scipy.org/doc/ in either the 0.7 or 0.8dev docs. I see splmake in my 0.7.0 installation, though (in the interpolate/interpolate.py file). Is the function splmake deprecated? I also notice that it is called (in the code) in the interp1d class init function, but it's impossible to specify the type (e.g., a clamped spline). Is there a reason why the kind argument is not passed to splmake? Thanks, Jason -- Jason Grout

2 5

Crash with Report from scikits.timeseries
by Sergi Pons Freixes 02 Mar '10

02 Mar '10

Hi All, When creating a simple report, I get an error which I'm not able to find out the cause: #Load the modules: In [1]: import scikits.timeseries as ts In [2]: import scikits.timeseries.lib.reportlib as rl #Load data: In [3]: labels = ["Temp", "S", "O2(%)", "dinophyceae", "diatoms", "chlorophyceae", "resta", "pseudonitzschia", "prymnesiophyceae(cocolitoforal)", "dinophysis", "alexandrium_minutum", "alexandrium_catenella", "dinophysis_caudata", "dinophysis_rotundata", "dinophysis_sacculus", "event"] In [4]: cols = [0] + range(3,19) In [5]: len(labels) Out[5]: 16 In [6]: len(cols) Out[6]: 17 In [7]: tmp = ts.tsfromtxt("AS.txt", freq='D', delimiter='\t', skip_header=1, usecols=cols, datecols=0, missing_values='nan' ) In [8]: tmp.data.shape Out[8]: (711,) #Create Report In [9]: report = rl.Report(tmp) In [10]: report() --------------------------------------------------------------------------- IndexError Traceback (most recent call last) /home/sergi/Dropbox/doctorat/dades_alfacs/<ipython console> in <module>() /usr/local/lib/python2.6/dist-packages/scikits.timeseries-0.91.3-py2.6-linux-i686.egg/scikits/timeseries/lib/reportlib.pyc in __call__(self, *tseries, **kwargs) 344 [datefmt_func(d)] + \ 345 [fmt_func[i](ser.series[d - _sd]) \ --> 346 for i, ser in enumerate(tseries)] 347 ) 348 /usr/lib/python2.6/dist-packages/numpy/ma/core.pyc in __getitem__(self, indx) 2576 # msg = "Masked arrays must be filled before they can be used as indices!" 2577 # raise IndexError, msg -> 2578 dout = ndarray.__getitem__(ndarray.view(self, ndarray), indx) 2579 # We could directly use ndarray.__getitem__ on self... 2580 # But then we would have to modify __array_finalize__ to prevent the IndexError: index out of bounds Is there something wrong with my data (even it is loaded with no errors). Regards, Sergi -- Sergi Pons Freixes Ph.D. Student Marine Technology Unit Centre Mediterrani d'Investigacions Marines i Ambientals (CMIMA-CSIC) Pg. Marítim Barceloneta, 37-49 E-08003 Barcelona (Spain) Ph. +34 93 230 95 00 (ext. 1510) Fax. +34 93 230 95 55 spons(a)utm.csic.es http://www.utm.csic.es

1 0

Re: [SciPy-User] scikits.learn installation
by David Trethewey 02 Mar '10

02 Mar '10

Problem is sorted now, thanks for your help.

1 0