Mailman 3 September 2013 - NumPy-Discussion

Re: [Numpy-discussion] Valid algorithm for generating a 3D Wiener Process?
by David Goldsmith 27 Sep '13

27 Sep '13

Thanks, guys. Yeah, I realized the problem w/ the uniform-increment-variable-direction approach this morning: physically, it ignores the fact that the particles hitting the particle being tracked are going to have a distribution of momentum, not all the same, just varying in direction. But I don't quite understand Warren's observation: "the 'angles' that describe the position undergo a random walk [actually, it would seem that they don't, since they too fail the varying-as-white-noise test], so the particle tends to move in the same direction over short intervals"--is this just another way of saying that, since I was varying the angles by -1, 0, or 1 unit each time, the simulation is susceptible to "unnaturally" long strings of -1, 0, or 1 increments? Thanks again, DG

4 4

Re: [Numpy-discussion] invalid correlation coefficient from np.ma.corrcoef
by Faraz Mirzaei 26 Sep '13

26 Sep '13

Thanks Josef and Nathaniel for your responses. In the application that I have, I don't use the correlation coefficient matrix as a whole (so I don't care if it is PSD or not). I simply read the off-diagonal elements for pair-wise correlation coefficients. I use the pairwise correlation coefficient to test if the data from various sources (i.e., rows of the matrix), agree with each other when present. Right now, I use, ma.corrcoef( x[ i, :] , x[ j, :]) and read the off-diagonal element in a loop over i and j. It is just a bit uglier than calling ma.corrcoef(x). At least for my application, truncation to -1 or +1 (or scaling such that largest values becomes 1 etc) is completely wrong, since it would imply that the two sources completely agree with each other (factoring out a minus sign), which may not the case. For example, consider the first and last rows of the example I provided: >>> print x_ma [[ 7 -4 -1 -- -3 -2] [ 6 -3 0 4 0 5] [-4 -- 7 5 -- --] [-- 5 -- 0 1 4]] >>> np.ma.corrcoef(x_ma)[0,3] -1.6813456149534147 On the other hand, if we supply only the first and third row to the function, we get: >>> np.ma.corrcoef(x_ma[0,:], x_ma[3,:]) masked_array(data = [[1.0 -0.240192230708] [-0.240192230708 1.0]], mask = [[False False] [False False]], fill_value = 1e+20) Interestingly, this is the same as what pandas results as the [3,0] element of the correlation coefficient matrix, and it is equal to pair-wise deletion result: >>> np.corrcoef([-4, -3, -2], [5, 1, 4]) #Note that this is NOT np.ma.corrcoef >>> array([[ 1. , -0.24019223], [-0.24019223, 1. ]]) Also, I don't know why the ma.corrcoef results Josef has mentioned are different than mine. In particular, Josef reports element [2, 0] of the ma.corrcoef result to be -1.19, but I get -- (i.e., missing and masked, probably due to too few samples available). Josef: are you sure that you have entered the example values correctly into python? Along the same lines, the results that Nathaniel has posted from R are different, since the input is not a masked matrix I guess (please note that in the original example, I had masked values less than or equal to -5). In any case, I think the correlation coefficient between two rows of a matrix should not depend on what other rows are supplied. In other words, np.ma.corrcoef(x_ma)[0,3] should be equal to np.ma.corrcoef(x_ma[0,:], x_ma[3,:])[0,1] (which apparently happens to be what pandas reports). This change would need recomputing the mean for every pair-wise coefficient calculation, but since we are computing cross products O(n^2) times, the overall big-O complexity won't change. And please don't remove this functionality. I will volunteer to fix it however we decide :) We can just clarify the behavior in the documentation. Thanks, Faraz

3 2

At.: use less RAM memory and increase execution speed
by Josè Luis Mietta 26 Sep '13

26 Sep '13

Hi experts! I wanna use less RAM memory in my Monte Carlo simulations. In my algorithm I use numpy arrays and xrange() function. I hear that I can reduce RAM used in my lagorithm if I do the next: 1) replace xrange() for range(). 2) replace numpya arrays for python lists 3) use reset() function for deleting useless arrays. Is that true? In adition, I wanna increase execution speed of my code (I use numpy and SciPy functions). How can I apply Cython? Will it help? Please help. Thanks a lot!!

3 2

invalid correlation coefficient from np.ma.corrcoef
by Faraz Mirzaei 26 Sep '13

26 Sep '13

Hi everyone, I'm using np.ma.corrcoef to compute the correlation coefficients among rows of a masked matrix, where the masked elements are the missing data. I've observed that in some cases, the np.ma.corrcoef gives invalid coefficients that are greater than 1 or less than -1. Here's an example: x = array([[ 7, -4, -1, -7, -3, -2], [ 6, -3, 0, 4, 0, 5], [-4, -5, 7, 5, -7, -7], [-5, 5, -8, 0, 1, 4]]) x_ma = np.ma.masked_less_equal(x , -5) C = np.round(np.ma.corrcoef(x_ma), 2) print C [[1.0 0.73 -- -1.68] [0.73 1.0 -0.86 -0.38] [-- -0.86 1.0 --] [-1.68 -0.38 -- 1.0]] As you can see, the [0,3] element is -1.68 which is not a valid correlation coefficient. (Valid correlation coefficients should be between -1 and 1). I looked at the code for np.ma.corrcoef, and this behavior seems to be due to the way that mean values of the rows of the input matrix are computed and subtracted from them. Apparently, the mean value is individually computed for each row, without masking the elements corresponding to the masked elements of the other row of the matrix, with respect to which the correlation coefficient is being computed. I guess the right way should be to recompute the mean value for each row every time that a correlation coefficient is being computed for two rows after propagating pair-wise masked values. Please let me know what you think. Thanks, Faraz

3 8

Valid algorithm for generating a 3D Wiener Process?
by David Goldsmith 25 Sep '13

25 Sep '13

Is this a valid algorithm for generating a 3D Wiener process? (When I graph the results, they certainly look like potential Brownian motion tracks.) def Wiener3D(incr, N): r = incr*(R.randint(3, size=(N,))-1) r[0] = 0 r = r.cumsum() t = 2*np.pi*incr*(R.randint(3, size=(N,))-1) t[0] = 0 t = t.cumsum() p = np.pi*incr*(R.randint(3, size=(N,))-1) p[0] = 0 p = p.cumsum() x = r*np.cos(t)*np.sin(p) y = r*np.sin(t)*np.sin(p) z = r*np.cos(p) return np.array((x,y,z)).T Thanks! DG

3 3

Heads up re: new memory tracking infrastructure in python 3.4
by Nathaniel Smith 25 Sep '13

25 Sep '13

Hi all, Possibly of interest to some folks here, it looks like in Python 3.4 there'll be a generic system for hooking and tracking memory allocations: http://www.python.org/dev/peps/pep-0445/ I'm not planning to do anything with this information myself, but if anyone's been thinking about numpy memory profiling, the new numpy malloc hooks, etc., then it might make sense to see how we could integrate with this more general system. -n

1 0

Numpy newbie question: matrix creation
by Edmondo Porcu 25 Sep '13

25 Sep '13

Dear all, I am a Newbie with Numpy and I would need some advice on how to create a matrix with certain characteristics : - Every entry should be minimum 0 maximum 1 with a step of 0.1 (legal values are 0,0.1,0.2,0.3 etc) - The number of columns of the matrix is a parameter of this matrix creation algorithm - Only the rows where the sum is 1 must be kept Would great appreciate your advice and suggestions Best Regards Edmondo

4 5

Removal of numarray and oldnumeric packages.
by Charles R Harris 24 Sep '13

24 Sep '13

Hi All, I have gotten no feedback on the removal of the numarray and oldnumeric packages. Consequently the removal will take place on 9/28. Scream now or never... Chuck

5 9

Pull request review #3770: Trapezoidal distribution
by Jeremy Hetzel 24 Sep '13

24 Sep '13

I've added a trapezoidal distribution to numpy.random for consideration, pull request 3770: https://github.com/numpy/numpy/pull/3770 Similar to the triangular distribution, the trapezoidal distribution may be used where the underlying distribution is not known, but some knowledge of the limits and mode exists. The trapezoidal distribution generalizes the triangular distribution by allowing the modal values to be expressed as a range instead of a point estimate. The trapezoidal distribution implemented, known as the "generalized trapezoidal distribution," has three additional parameters: growth, decay, and boundary ratio. Adjusting these from the default values create trapezoidal-like distributions with non-linear behavior. Examples can be seen in an R vignette ( http://cran.r-project.org/web/packages/trapezoid/vignettes/trapezoid.pdf ), as well as these papers by J.R. van Dorp and colleagues: 1) van Dorp, J. R. and Kotz, S. (2003) Generalized trapezoidal distributions. Metrika. 58(1):85–97. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/Metrika2003VanDo… 2) van Dorp, J. R., Rambaud, S.C., Perez, J. G., and Pleguezuelo, R. H. (2007) An elicitation procedure for the generalized trapezoidal distribution with a uniform central stage. Decision Analysis Journal. 4:156–166. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/DA2007.pdf The docstring for the proposed numpy.random.trapezoidal() is as follows: """ trapezoidal(left, mode1, mode2, right, size=None, m=2, n=2, alpha=1) Draw samples from the generalized trapezoidal distribution. The trapezoidal distribution is defined by minimum (``left``), lower mode (``mode1``), upper mode (``mode1``), and maximum (``right``) parameters. The generalized trapezoidal distribution adds three more parameters: the growth rate (``m``), decay rate (``n``), and boundary ratio (``alpha``) parameters. The generalized trapezoidal distribution simplifies to the trapezoidal distribution when ``m = n = 2`` and ``alpha = 1``. It further simplifies to a triangular distribution when ``mode1 == mode2``. Parameters ---------- left : scalar Lower limit. mode1 : scalar The value where the first peak of the distribution occurs. The value should fulfill the condition ``left <= mode1 <= mode2``. mode2 : scalar The value where the first peak of the distribution occurs. The value should fulfill the condition ``mode1 <= mode2 <= right``. right : scalar Upper limit, should be larger than or equal to `mode2`. size : int or tuple of ints, optional Output shape. Default is None, in which case a single value is returned. m : scalar, optional Growth parameter. n : scalar, optional Decay parameter. alpha : scalar, optional Boundary ratio parameter. Returns ------- samples : ndarray or scalar The returned samples all lie in the interval [left, right]. Notes ----- With ``left``, ``mode1``, ``mode2``, ``right``, ``m``, ``n``, and ``alpha`` parametrized as :math:`a, b, c, d, m, n, \\text{ and } \\alpha`, respectively, the probability density function for the generalized trapezoidal distribution is .. math:: f{\\scriptscriptstyle X}(x\mid\theta) = \\mathcal{C}(\\Theta) \\times \\begin{cases} \\alpha \\left(\\frac{x - \\alpha}{b - \\alpha} \\right)^{m - 1}, & \\text{for } a \\leq x < b \\\\ (1 - \\alpha) \\left(\frac{x - b}{c - b} \\right) + \\alpha, & \\text{for } b \\leq x < c \\\\ \\left(\\frac{d - x}{d - c} \\right)^{n-1}, & \\text{for } c \\leq x \\leq d \\end{cases} with the normalizing constant :math:`\\mathcal{C}(\\Theta)` defined as ..math:: \\mathcal{C}(\\Theta) = \\frac{2mn} {2 \\alpha \\left(b - a\\right) n + \\left(\\alpha + 1 \\right) \\left(c - b \\right)mn + 2 \\left(d - c \\right)m} and where the parameter vector :math:`\\Theta = \\{a, b, c, d, m, n, \\alpha \\}, \\text{ } a \\leq b \\leq c \\leq d, \\text{ and } m, n, \\alpha >0`. Similar to the triangular distribution, the trapezoidal distribution may be used where the underlying distribution is not known, but some knowledge of the limits and mode exists. The trapezoidal distribution generalizes the triangular distribution by allowing the modal values to be expressed as a range instead of a point estimate. The growth, decay, and boundary ratio parameters of the generalized trapezoidal distribution further allow for non-linear behavior to be specified. References ---------- .. [1] van Dorp, J. R. and Kotz, S. (2003) Generalized trapezoidal distributions. Metrika. 58(1):85–97. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/Metrika2003VanDo… .. [2] van Dorp, J. R., Rambaud, S.C., Perez, J. G., and Pleguezuelo, R. H. (2007) An elicitation proce-dure for the generalized trapezoidal distribution with a uniform central stage. Decision AnalysisJournal. 4:156–166. Preprint available: http://www.seas.gwu.edu/~dorpjr/Publications/JournalPapers/DA2007.pdf Examples -------- Draw values from the distribution and plot the histogram: >>> import matplotlib.pyplot as plt >>> h = plt.hist(np.random.triangular(0, 0.25, 0.75, 1, 100000), bins=200, ... normed=True) >>> plt.show() """ I am unsure if NumPy encourages incorporation of new distributions into numpy.random or instead into separate modules, but found the exercise to be helpful regardless. Thanks, Jeremy

3 4

Binary releases
by Charles R Harris 24 Sep '13

24 Sep '13

Hi All, Numpy 1.8 is about ready for an rc1, which brings up the question of which binary builds so put up on sourceforge. For Windows maybe 1. 32 bit windows, python 2.6, 2.7, 3.2, 3.3, compiled with MSVC 2. 64 bit windows, python 2.6, 2.7, 3.2, 3.3, compiled with MSVC, linked with MKL I think these should be good for both windows 7 and window 8, correct me if I am wrong. For Mac there is first the question of OS X versions, (10.5?), 10.6, 10.7, 10.8. I don't know if some builds work on more than one OS X version. The 10.5 version is a bit harder to come by than 10.6 and up. It looks like 10.9 is coming up, but it isn't out yet. I have no idea what Python version to match up these, but assuming all of them, then 1. OS X 10.6 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. 2. OS X 10.7 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. 3. OS X 10.8 python 2.6, 2.7, 3.2, 3.3, compiled with native compiler, linked with Accelerate. That seems like a lot. It is fairly easy to compile from source on the mac these days, are all those binary packages really needed? I don't know what I am doing with the binary stuff, so any suggestions are welcome. For building it is possible to set up a Macintosh with vmware fusion to handle all of these, but there is time and money involved in that. Any one who is already capable of doing some of these builds is welcome to volunteer. Note, Intel has offered MKL licenses to the open source projects under the NumFocus umbrella, but I don't know of anyone who has taken advantage of that yet or how much time it will take to go through the undoubtedly needed paper work, but I would like to get some of those licenses and move to MSVC so that we can stop rolling around gasping in pain whenever it comes to window builds. Then there is Fortran for windows... Thoughts? Chuck

11 31