<p dir="auto">My issue is about a funny behavior of scipy function integrate.quad.</p>
<p dir="auto">Integrating a gaussian function I get wrong values depending on the integration interval. See below.<br>
in [3] I compute the integral on [-100, 0.5] and I get the expected value.</p>
<p dir="auto">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.</p>
<p dir="auto">In [5] and [6] I show that something happens between -997 and -998</p>
<p dir="auto">In [7] I extend further the integration interval and the returned value is still wrong (I've tried several values with the same result.</p>
<p dir="auto">In [8] I show that the integral is correct if I integrate from -infinity</p>
<h4 dir="auto">Reproducing code example:</h4>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="$ 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)
"><pre class="notranslate"><code class="notranslate">$ 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)

</code></pre></div>
<h4 dir="auto">Error message:</h4>
<p dir="auto">No error message. Wrong result</p>
<h4 dir="auto">Scipy/Numpy/Python version information:</h4>
<p dir="auto">Hi,<br>
I was just exploring a difference I noticed between numpy's std() function and scipy.stats' nanstd() function.</p>
<p dir="auto">numpy.std is called like this:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False)"><pre class="notranslate"><code class="notranslate">numpy.std(a, axis=None, dtype=None, out=None, ddof=0, keepdims=False)
</code></pre></div>
<p dir="auto">and scipy.stats.nanstd is called like this:</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="nanstd(x, axis=0, bias=False)"><pre class="notranslate"><code class="notranslate">nanstd(x, axis=0, bias=False)
</code></pre></div>
<p dir="auto">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.</p>
<p dir="auto">This could be fixed like so, however I can't think of a clever way to keep compatibility with the old bias keyword. Opinions?</p>
<div class="snippet-clipboard-content notranslate position-relative overflow-auto" data-snippet-clipboard-copy-content="def nanstd(x, axis=0, ddof = 0):
    &quot;&quot;&quot;
    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
    --------
    &gt;&gt;&gt; from scipy import stats
    &gt;&gt;&gt; a = np.arange(10, dtype=float)
    &gt;&gt;&gt; a[1:3] = np.nan
    &gt;&gt;&gt; np.std(a)
    nan
    &gt;&gt;&gt; stats.nanstd(a)
    2.9154759474226504
    &gt;&gt;&gt; stats.nanstd(a.reshape(2, 5), axis=1)
    array([ 2.0817, 1.5811])
    &gt;&gt;&gt; stats.nanstd(a.reshape(2, 5), axis=None)
    2.9154759474226504

&quot;&quot;&quot;
    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)"><pre class="notranslate"><code class="notranslate">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
    --------
    &gt;&gt;&gt; from scipy import stats
    &gt;&gt;&gt; a = np.arange(10, dtype=float)
    &gt;&gt;&gt; a[1:3] = np.nan
    &gt;&gt;&gt; np.std(a)
    nan
    &gt;&gt;&gt; stats.nanstd(a)
    2.9154759474226504
    &gt;&gt;&gt; stats.nanstd(a.reshape(2, 5), axis=1)
    array([ 2.0817, 1.5811])
    &gt;&gt;&gt; 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)
</code></pre></div>