Issue 1:
My issue is about a funny behavior of scipy function integrate.quad.

Integrating a gaussian function I get wrong values depending on the integration interval. See below.

in [3] I compute the integral on [-100, 0.5] and I get the expected value.

In [4] I compute the same integral on [-1000,0.5]. I expect the same value because there is no co ntribute to the integral sum extending the interval to the left. But I get a wrong value.

In [5] and [6] I show that something happens between -997 and -998

In [7] I extend further the integration interval and the returned value is still wrong (I've tried several values with the same result.

In [8] I show that the integral is correct if I integrate from -infinity

Reproducing code example:

$ ipython --pylab
Python 3.6.9 (default, Nov 7 2019, 10:44:02)
Type 'copyright', 'credits' or 'license' for more information
IPython 7.5.0 -- An enhanced Interactive Python. Type '?' for help.
In [1]: import scipy.integrate as integrate

In [2]: from scipy.stats import norm

In [3]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-100,0.5)
Out[3]: (0.22347533732644612, 3.1634786288170372e-12)

In [4]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-1000,0.5)
Out[4]: (4.75208186635514e-09, 9.44858985141507e-09)

In [5]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-987,0.5)
Out[5]: (0.22347533732644612, 1.8395174387908237e-10)

In [6]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-988,0.5)
Out[6]: (7.303423606392351e-09, 1.4521436311212504e-08)

In [7]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-2000,0.5)
Out[7]: (6.486552015239486e-32, 1.2897246147166321e-31)

In [8]: integrate.quad(lambda x:norm.pdf(x,0.64449662,0.19),-np.inf,0.5)
Out[8]: (0.2234753373264382, 1.598765108016799e-09)

Error message:

No error message. Wrong result

Scipy/Numpy/Python version information:

Issue 2:
Hi,

I was just exploring a difference I noticed between numpy's std() function and scipy.stats' nanstd() function.

numpy.std is called like this:

numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False)

and scipy.stats.nanstd is called like this:

nanstd(x, axis=0, bias=False)

The setting bias = False is equivalent to setting ddof=1 in numpy.std, leading to the calls numpy.std(x) and scipy.stats.nanstd(x) to return different values, especially for small arrays. I think these should by default return the same value.

This could be fixed like so, however I can't think of a clever way to keep compatibility with the old bias keyword. Opinions?

def nanstd(x, axis=0, ddof = 0):
"""
Compute the standard deviation over the given axis, ignoring nans.

Parameters
----------
x : array_like
Input array.
axis : int or None, optional
Axis along which the standard deviation is computed. Default is 0.
If None, compute over the whole array `x`.
ddof : int, optional
Means Delta Degrees of Freedom. The divisor used in calculations is
N - ddof, where N represents the number of elements.
By default ddof is zero.

Returns
-------
s : float
The standard deviation.

See Also
--------
nanmean, nanmedian

Examples
--------
>>> from scipy import stats
>>> a = np.arange(10, dtype=float)
>>> a[1:3] = np.nan
>>> np.std(a)
nan
>>> stats.nanstd(a)
2.9154759474226504
>>> stats.nanstd(a.reshape(2, 5), axis=1)
array([ 2.0817, 1.5811])
>>> stats.nanstd(a.reshape(2, 5), axis=None)
2.9154759474226504

"""
x, axis = _chk_asarray(x,axis)
x = x.copy()
Norig = x.shape[axis]

Nnan = np.sum(np.isnan(x),axis)*1.0
n = Norig - Nnan

x[np.isnan(x)] = 0.
m1 = np.sum(x,axis)/n

if axis:
d = (x - np.expand_dims(m1, axis))**2.0
else:
d = (x - m1)**2.0

m2 = np.sum(d,axis)-(m1*m1)*Nnan

m2c = m2 / (n - ddof)

return np.sqrt(m2c)