Mailman 3 September 2013 - SciPy-User

stats.gaussian_kde sensitive to outliers
by David Soukal Sept. 29, 2013

Sept. 29, 2013

Hello folks, I was studying kernel density estimation in "The Elements of Statistical Learning" and playing with stats.gaussian_kde() to confirm my understanding. I found that, contrary to my expectation, the estimates seem very sensitive to outliers. The code at the end of this email shows this visually. It generates n=1000 points from normal distribution (d) and adds an outlier to it (d_o). As you can see in the graph the estimated density with the outlier looks very different. Inspecting the estimated density at 0 clearly shows that a lot of mass has shifted when the outlier is present. print(kde(0)) print(kde_o(0)) [ 0.00389799] [ 0.0009715] I thought this could be because of the automatic bandwidth estimate but it doesn't seem so. I was surprised at this behavior. I thought that an outlier should have little or no impact on the density elsewhere since the Gaussian kernel is exponentially decaying. Fitting this in R did produce a more stable estimate. Can you, please, help me understand what I'm missing? Is gaussian_kde() doing something different than a standard smoothing where kernels are placed at the sample points and added up? THANKS, David # import matplotlib.pylab as plt import numpy as np import scipy.stats as stats # generate data with and without an outlier n = 1000 d = np.random.normal(0, 100, n) outlier = 50000 d_o = np.insert(d, 0, outlier) # estimate the densities kde = stats.gaussian_kde(d) kde_o = stats.gaussian_kde(d_o) # plot xs = np.linspace(-400, 400, 1000) fg, ax = plt.subplots(nrows = 1, ncols = 1, figsize = (15, 5)) ax.hist(d, bins = 100, normed = True, alpha = 0.5, label = 'Histogram') ax.plot(xs, kde(xs), linewidth = 2, alpha = 0.8, label = 'Gaussian KDE, no outlier') ax.plot(xs, kde_o(xs), linewidth = 2, alpha = 0.8, label = 'Gaussian KDE, outlier') ax.legend(loc = 'best') R code for reference. library(ggplot2) c1 <- rnorm(1000, 0, 100) d1 <- data.frame(x = c1) c2 <- c(c1, 50000) d2 <- data.frame(x = c2) kde1 <- density(d1$x, kernel = 'gaussian') kde2 <- density(d2$x, kernel = 'gaussian') ggplot() + geom_density(data = d1, aes(x = x), color = 'red') + geom_density(data = d2, aes(x = x), color = 'green') + coord_cartesian(xlim=c(-400, 400))

2 3

How do I cite SciPy and NumPy, when I used them for a paper?
by Yuxiang Wang Sept. 28, 2013

Sept. 28, 2013

Dear all, A quick question - how do I cite SciPy and NumPy, when I used them for a paper? Do I just say them without referencing, or should I refer to the www.scipy.org? Thanks in advance :) -Shawn -- Yuxiang "Shawn" Wang Gerling Research Lab University of Virginia yw5aj(a)virginia.edu +1 (434) 284-0836

3 3

Baffling error: ndarray += csc_matrix -> "ValueError: setting an array element with a sequence"
by Nathaniel Smith Sept. 27, 2013

Sept. 27, 2013

Hi all, I'm getting a very strange error, and unfortunately I can't seem to reproduce it *except* when running on Travis-CI, so my debugging tools are really limited. I'm wondering if anyone else has seen anything like this? The offending line of code is: File "/home/travis/virtualenv/python2.7_with_system_site_packages/local/lib/python2.7/site-packages/pyrerp/incremental_ls.py", line 323, in append_bottom_half self.xtx += xtx ValueError: setting an array element with a sequence. And debug prints reveal that in the case that causes the error, 'self.xtx' is an ndarray that prints as: [[ 0. 0.] [ 0. 0.]] while 'xtx' is a scipy.sparse.csc.csc_matrix that prints as: (1, 0) 45.0 (0, 0) 10.0 (1, 1) 285.0 (0, 1) 45.0 On my laptop (Ubuntu 12.10, 64-bit, numpy 1.7.1, scipy 0.12.0), the same test passes fine. Travis is Ubuntu 12.04, 32-bit, numpy 1.7.1, scipy 0.12.0, so I guess in *principle* it could be a 32-bit/64-bit thing, but that's just a wild guess. Any ideas? -n P.S.: if anyone's curious, the full build log is here: https://s3.amazonaws.com/archive.travis-ci.org/jobs/11735709/log.txt and the code is here: https://github.com/njsmith/pyrerp/blob/master/pyrerp/incremental_ls.py#L323

3 6

Minimization of a function that can return nan
by federico vaggi Sept. 27, 2013

Sept. 27, 2013

Hi everyone, I'm using fmin_l_bfgs_b to find the best fit parameters for a system of ordinary differential equations. Since the system has multiple local minima, I found that it's very convenient to divide the parameter space in hypercubes, and perform multiple local minimizations. Occasionally, for certain regions of the parameter space, the numerical integration routine fails. What's the correct way to handle this gracefully? Capture the exception, and return inf as the cost function value, and set the gradient to equal to all nan? Federico

1 0

Re: [SciPy-User] Sparse dot differs from numpy dot by 10^5
by Wagner Sebastian Sept. 27, 2013

Sept. 27, 2013

So, I can't get rid of this issue, because this error does not lie within scipy or numpy but in the design of floating point arithmetic? At the end of the calculation (after 3 more operations) only two promille out of 1470 values do not differ, all others do differ from the third digit on (median of the difference). Is there any chance to know which one of the two results is "more correct"? Regards, Sebastian -----Ursprüngliche Nachricht----- Von: scipy-user-bounces(a)scipy.org [mailto:scipy-user-bounces@scipy.org] Im Auftrag von Emanuele Olivetti Gesendet: Freitag, 27. September 2013 09:37 An: scipy-user(a)scipy.org Betreff: Re: [SciPy-User] Sparse dot differs from numpy dot by 10^5 On 09/27/2013 09:31 AM, Jerome Kieffer wrote: > On Fri, 27 Sep 2013 07:12:06 +0000 > Wagner Sebastian <Sebastian.Wagner.fl(a)ait.ac.at> wrote: > >> When changing some algorithms from dense matrices to sparse matrices I stumbled over a discrepancy between the sparse-version of dot and numpy's dot. All open issues mentioning 'dot' do not apply. > Floating point addition is not commutable ... the order you are > summing is important (in silico) while it is not on the paper. > If that is the problem, this is a detailed description: http://docs.oracle.com/cd/E19957-01/806-3568/ncg_goldberg.html Best, Emanuele _______________________________________________ SciPy-User mailing list SciPy-User(a)scipy.org http://mail.scipy.org/mailman/listinfo/scipy-user

2 1

Sparse dot differs from numpy dot by 10^5
by Wagner Sebastian Sept. 27, 2013

Sept. 27, 2013

Dear Scipy-Users, When changing some algorithms from dense matrices to sparse matrices I stumbled over a discrepancy between the sparse-version of dot and numpy's dot. All open issues mentioning 'dot' do not apply. The data-source is a square matrix with N=1470, which is originally created by the constructor of scipy.sparse.csc_matrix from a dense matrix (You can recreate it by converting the matrix first to dense and then back to csc). I compute the dot-product of the matrix with itself transposed and some of the dot products (28 exactly) are different from the np.dot by a factor of 10^5. I saved the csc-matrix to a npz-file and attached it. A code sample to reproduce the effect is here: sparsecsc = np.load('sparsecsc.npz') K = sparsecsc['K'][()] K[249].dot(K[:,251]).A # gives -9.61216512e+08 np.dot(K[249].A,K[:,251].A # gives -9.61150976e+08 I located the diverging cells by using the equations behind np.allclose(): def close(a, b, rtol=1e-05, atol=1e-08): c = np.absolute(a-b) <= (atol + rtol * np.absolute(b)) d = np.absolute(b-a) <= (atol + rtol * np.absolute(a)) return c | d c = close(K.dot(K.T).A, np.dot(K.A, K.T.A)) c = ~c c.nonzero() This gives pairs of 28 indices which diverge. Or to get all diverging cells: K.dot(K.T).A[c] np.dot(K.A, K.T.A)[c] # or the derived differences: K.dot(K.T).A[c] - np.dot(K.A, K.T.A)[c] # average of them: 541257 Does anyone have an idea what's going on here? Regards, Sebastian

3 2

Re: [SciPy-User] [cython-users] scipy cblas return value
by Radim Řehůřek Sept. 24, 2013

Sept. 24, 2013

It's actually more serious than I though. This "double instead of float" apparently fails on some systems while manifesting itself as a weird slow down: https://groups.google.com/d/msg/gensim/6B99kdtHOZo/FSss6GgzE5IJ We've confirmed with the person who reported it that the problem is only in sdot. Changing the sdot return value to `float` fixes it. SciPy version doesn't seem to affect this (tried both 0.12 and 0.13b on the failing machine). CC to Scipy people: the code is using a protected pointer from scipy.linalg.blas.cblas, like so: https://gist.github.com/anonymous/6659007 Why is the required return value sometimes float and sometimes double? Is there a better/safer way to achieve calling numpy/scipy BLAS from cython? Cheers, Radim On Tuesday, September 24, 2013 6:58:20 AM UTC+2, Robert Bradshaw wrote: > > _cpointer seems to be woefully undocumented, but it looks like it is a > C function actually returning a double not a float. > > On Sun, Sep 22, 2013 at 4:54 AM, Radim Řehůřek <m...(a)radimrehurek.com<javascript:>> > wrote: > > Hello all, > > > > I'm trying to call sdot from scipy.linalg.blas.cblas from Cython, but it > > only works if I declare the return type as double (not float): > > > > https://gist.github.com/anonymous/6659007 > > > > I don't know if this is a Cython issue or a SciPy issue or a BLAS > (vecLib) > > issue or a chair-keyboard issue. Could someone more knowledgeable > elucidate > > please? > > > > Cheers, > > Radim, Cython 0.19.1, SciPy 0.12.0, Python 2.7.1, NumPy & SciPy BLAS > from OS > > X Accelerate > > > > > > -- > > > > --- > > You received this message because you are subscribed to the Google > Groups > > "cython-users" group. > > To unsubscribe from this group and stop receiving emails from it, send > an > > email to cython-users...(a)googlegroups.com <javascript:>. > > For more options, visit https://groups.google.com/groups/opt_out. >

2 1

Akima interpolation?
by Andreas Hilboll Sept. 22, 2013

Sept. 22, 2013

I'm looking for 1d Akima interpolation, and couldn't find it in Scipy. 1.) Did I miss something? 2.) Is there interest in having Akima interpolation in Scipy? 3.) Does anyone know of a good, fast interpolation I should investigate for inclusion in Scipy? Cheers, Andreas. -- Andreas Hilboll PhD Student Institute of Environmental Physics University of Bremen U3145 Otto-Hahn-Allee 1 D-28359 Bremen Germany +49(0)421 218 62133 (phone) +49(0)421 218 98 62133 (fax) http://www.iup.uni-bremen.de/~hilboll

2 1

discrete integration code ?
by josef.pktd＠gmail.com Sept. 22, 2013

Sept. 22, 2013

general question: Is there something like integrate.quad for discrete integration sum of f(i) over integers (with or without finite bounds) the main point is to figure out some termination condition scipy.stats has a "home made" function, but we ran into a problem where the summation stops too soon. Are there any good rules when to stop summing terms? usually we just benefit from the `special` functions. our special case f(i) = g(i) * p(i) where p(i) is the probability mass function, and g(i) could be any function Josef

1 0

Small error in example of documentation of signal.freqz?
by Alejandro Weinstein Sept. 22, 2013

Sept. 22, 2013

Hi: The documentation of signal.freqz [1] has an example that plots the frequency response of a filter. The y axis is plotted in a logaritmic scale and labeled as "decibels": plt.semilogy(w, np.abs(h), 'b') plt.ylabel('Amplitude (dB)', color='b') I think that to be in dB, the plot should be in a linear scale and the magnitude of the response computed as 20*np.log10(np.abs(h)): plt.plot(w, 20*np.log10(np.abs(h)), 'b') plt.ylabel('Amplitude (dB)', color='b') Any comments? Alejandro [1] http://docs.scipy.org/doc/scipy/reference/generated/scipy.signal.freqz.html

2 1