Repository: dask/dask
Pull request: Parallel variance computation for dataframes
Fixes #4233

- [x] Tests added / passed
- [x] Passes `flake8 dask`

Re-using array variance computation logic in dataframe.

Problem statement:
Dask Series std function returns NaN for valid array
When executing `memory_usage_series.std().compute()` where `memory_usage_series = task_dataframe["memory_consumption"]` I got a Warning.

```
/home/lfdversluis/miniconda3/envs/format2/lib/python3.6/site-packages/dask/compatibility.py:95: RuntimeWarning: invalid value encountered in sqrt
  return func(*args)
```

and the output is NaN.

I checked the Series and all the values are fine. They are np.int64 integers and there are no `NaN`, or `None` values in the array. If I now parse the same array directly using numpy from a pandas dataframe (on the same data, working on a subset right now), `print(task_df["memory_consumption"].values.std(dtype=np.float64))` I do get the correct output. I cannot find any reasons in the document that could case this bug. I also tried passing `dtype=np.float64` in the Dask std call, but I get the same warning and output.

Hints from issue discussion:
Thanks for the bug report @lfdversluis .  If you're able to provide a minimal reproducible example that would help maintainers be able to identify what is wrong more quickly.
Also possibly a repeat of https://github.com/dask/dask/issues/3906
Sure thing. I don't think it is related to #3906, there are no gaps AFAIK.
Here is a simple example:

```
import pandas as pd
import numpy as np
import dask.dataframe as dd

pandas_df = pd.read_csv("test.csv", low_memory=False, memory_map=True, engine='c', dtype={"memory_consumption": np.int64}, header=0,
                delimiter=",")

print(pandas_df['memory_consumption'].std())

dask_df = dd.from_pandas(pandas_df, npartitions=1)
print(dask_df['memory_consumption'].std().compute())
```

`test.csv`:  [test.zip](https://github.com/dask/dask/files/2607515/test.zip)
I could reproduce it as well.

This code returns NaN:
```python 
ddf = dd.from_pandas(pd.DataFrame({'a': [13094304034, 3489385935, 100006774]}), npartitions=2)
ddf['a'].std().compute()
```

Warning:
`compatibility.py:95: RuntimeWarning: invalid value encountered in sqrt 
return func(*args)`

When I debugged there was a negative variance passed to np.sqrt function.

Looks like it overflows in var function:
```python 
x2 = 1.0 * (num ** 2).sum(skipna=skipna, split_every=split_every)
```
Maybe smth like this can be used: https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
@jcrist did something similar for Dask Array a while ago in
https://github.com/dask/dask/pull/297 and it seems to have worked out
well.  Doing the same for Dask DataFrame would be a nice (and possibly fun)
contribution.

On Tue, May 7, 2019 at 4:38 PM Ksenia Bobrova <notifications@github.com>
wrote:

> I could repro it as well.
>
> This code returns NaN:
>
> ddf = dd.from_pandas(pd.DataFrame({'a': [13094304034, 3489385935, 100006774]}), npartitions=2)
> ddf['a'].std().compute()
>
> Warning:
> compatibility.py:95: RuntimeWarning: invalid value encountered in sqrt
> return func(*args)
>
> When I debugged there was a negative variance passed to np.sqrt function.
>
> Looks like it overflows in var function:
>
> x2 = 1.0 * (num ** 2).sum(skipna=skipna, split_every=split_every)
>
> Maybe smth like this can be used:
> https://en.wikipedia.org/wiki/Algorithms_for_calculating_variance#Parallel_algorithm
>
> —
> You are receiving this because you commented.
> Reply to this email directly, view it on GitHub
> <https://github.com/dask/dask/issues/4233#issuecomment-490264837>, or mute
> the thread
> <https://github.com/notifications/unsubscribe-auth/AACKZTHQWIGXOA472HNAOA3PUHZFPANCNFSM4GFVV4CA>
> .
>

@mrocklin for dask array it works, indeed :) thank you, looks interesting.

Commit message: Adding support of timedelta64

Patch to review:
diff --git a/dask/array/reductions.py b/dask/array/reductions.py
index e22489a1d3b..f093ac4249b 100644
--- a/dask/array/reductions.py
+++ b/dask/array/reductions.py
@@ -6,7 +6,7 @@
 from math import factorial, log, ceil
 
 import numpy as np
-from numbers import Integral
+from numbers import Integral, Number
 
 from toolz import compose, partition_all, get, accumulate, pluck
 
@@ -457,7 +457,17 @@ def moment_agg(pairs, order=2, ddof=0, dtype='f8', sum=np.sum, axis=None, **kwar
     inner_term = divide(totals, ns, dtype=dtype) - mu
 
     M = _moment_helper(Ms, ns, inner_term, order, sum, axis, kwargs)
-    return divide(M, n.sum(axis=axis, **kwargs) - ddof, dtype=dtype)
+
+    denominator = n.sum(axis=axis, **kwargs) - ddof
+
+    # taking care of the edge case with empty or all-nans array with ddof > 0
+    if isinstance(denominator, Number):
+        if denominator < 0:
+            denominator = np.nan
+    else:
+        denominator[denominator < 0] = np.nan
+
+    return divide(M, denominator, dtype=dtype)
 
 
 def moment(a, order, axis=None, dtype=None, keepdims=False, ddof=0,
diff --git a/dask/dataframe/core.py b/dask/dataframe/core.py
index 793476e8954..02f8641a0a2 100644
--- a/dask/dataframe/core.py
+++ b/dask/dataframe/core.py
@@ -1461,17 +1461,76 @@ def var(self, axis=None, skipna=True, ddof=1, split_every=False, dtype=None, out
                                     axis=axis, skipna=skipna, ddof=ddof)
             return handle_out(out, result)
         else:
-            num = self._get_numeric_data()
-            x = 1.0 * num.sum(skipna=skipna, split_every=split_every)
-            x2 = 1.0 * (num ** 2).sum(skipna=skipna, split_every=split_every)
-            n = num.count(split_every=split_every)
-            name = self._token_prefix + 'var'
-            result = map_partitions(methods.var_aggregate, x2, x, n,
-                                    token=name, meta=meta, ddof=ddof)
+            if self.ndim == 1:
+                result = self._var_1d(self, skipna, ddof, split_every)
+            elif len(self._meta_nonempty.select_dtypes(include=[np.timedelta64]).columns) > 0:
+                result = self._var_mixed(skipna, ddof, split_every)
+            else:
+                result = self._var_numeric(skipna, ddof, split_every)
+
             if isinstance(self, DataFrame):
                 result.divisions = (min(self.columns), max(self.columns))
             return handle_out(out, result)
 
+    def _var_numeric(self, skipna=True, ddof=1, split_every=False):
+        num = self.select_dtypes(include=['number', 'bool'])
+
+        values_dtype = num.values.dtype
+        array_values = num.values
+
+        if not np.issubdtype(values_dtype, np.number):
+            array_values = num.values.astype('f8')
+
+        var = da.nanvar if skipna or skipna is None else da.var
+        array_var = var(array_values, axis=0, ddof=ddof, split_every=split_every)
+
+        name = self._token_prefix + 'var-numeric' + tokenize(num, split_every)
+        cols = num._meta.columns if is_dataframe_like(num) else None
+
+        var_shape = num._meta_nonempty.values.var(axis=0).shape
+        array_var_name = (array_var._name,) + (0,) * len(var_shape)
+
+        layer = {(name, 0): (methods.wrap_var_reduction, array_var_name, cols)}
+        graph = HighLevelGraph.from_collections(name, layer, dependencies=[array_var])
+
+        return new_dd_object(graph, name, num._meta_nonempty.var(), divisions=[None, None])
+
+    def _var_mixed(self, skipna=True, ddof=1, split_every=False):
+        data = self.select_dtypes(include=['number', np.timedelta64, 'bool'])
+
+        vars = [self._var_1d(data[col_idx], skipna, ddof, split_every) for col_idx in data._meta.columns]
+        var_names = [(v._name, 0) for v in vars]
+
+        name = self._token_prefix + 'var-mixed-' + tokenize(data, split_every)
+
+        layer = {(name, 0): (methods.wrap_var_reduction, var_names, data._meta.columns)}
+        graph = HighLevelGraph.from_collections(name, layer, dependencies=vars)
+
+        return new_dd_object(graph, name, data._meta_nonempty.var(), divisions=[None, None])
+
+    def _var_1d(self, column, skipna=True, ddof=1, split_every=False):
+        is_timedelta = is_timedelta64_dtype(column._meta)
+
+        if is_timedelta:
+            column = column.dropna().astype('i8')
+
+        if PANDAS_VERSION >= '0.24.0':
+            if pd.Int64Dtype.is_dtype(column._meta_nonempty):
+                column = column.astype('f8')
+
+        if not np.issubdtype(column.dtype, np.number):
+            column = column.astype('f8')
+
+        name = self._token_prefix + 'var-1d-' + tokenize(column, split_every)
+
+        var = da.nanvar if skipna or skipna is None else da.var
+        array_var = var(column.values, axis=0, ddof=ddof, split_every=split_every)
+
+        layer = {(name, 0): (methods.wrap_var_reduction, (array_var._name,), None)}
+        graph = HighLevelGraph.from_collections(name, layer, dependencies=[array_var])
+
+        return new_dd_object(graph, name, column._meta_nonempty.var(), divisions=[None, None])
+
     @derived_from(pd.DataFrame)
     def std(self, axis=None, skipna=True, ddof=1, split_every=False, dtype=None, out=None):
         axis = self._validate_axis(axis)
@@ -1503,6 +1562,7 @@ def sem(self, axis=None, skipna=None, ddof=1, split_every=False):
             n = num.count(split_every=split_every)
             name = self._token_prefix + 'sem'
             result = map_partitions(np.sqrt, v / n, meta=meta, token=name)
+
             if isinstance(self, DataFrame):
                 result.divisions = (min(self.columns), max(self.columns))
             return result
diff --git a/dask/dataframe/methods.py b/dask/dataframe/methods.py
index c31c52237f6..157e0106e37 100644
--- a/dask/dataframe/methods.py
+++ b/dask/dataframe/methods.py
@@ -121,16 +121,11 @@ def mean_aggregate(s, n):
         return np.float64(np.nan)
 
 
-def var_aggregate(x2, x, n, ddof):
-    try:
-        with warnings.catch_warnings(record=True):
-            warnings.simplefilter('always')
-            result = (x2 / n) - (x / n)**2
-        if ddof != 0:
-            result = result * n / (n - ddof)
-        return result
-    except ZeroDivisionError:
-        return np.float64(np.nan)
+def wrap_var_reduction(array_var, index):
+    if isinstance(array_var, np.ndarray) or isinstance(array_var, list):
+        return pd.Series(array_var, index=index)
+
+    return array_var
 
 
 def describe_aggregate(values):
diff --git a/dask/dataframe/tests/test_arithmetics_reduction.py b/dask/dataframe/tests/test_arithmetics_reduction.py
index d7a79e3e430..7fdf0ee42ca 100644
--- a/dask/dataframe/tests/test_arithmetics_reduction.py
+++ b/dask/dataframe/tests/test_arithmetics_reduction.py
@@ -607,13 +607,19 @@ def test_frame_series_arithmetic_methods():
 
 @pytest.mark.parametrize('split_every', [False, 2])
 def test_reductions(split_every):
-    dsk = {('x', 0): pd.DataFrame({'a': [1, 2, 3], 'b': [4, 5, 6]},
+    dsk = {('x', 0): pd.DataFrame({'a': [1, 2, 3],
+                                   'b': [4, 5, 6],
+                                   'c': [True, True, False]},
                                   index=[0, 1, 3]),
-           ('x', 1): pd.DataFrame({'a': [4, 5, 6], 'b': [3, 2, 1]},
+           ('x', 1): pd.DataFrame({'a': [4, 5, 6],
+                                   'b': [3, 2, 1],
+                                   'c': [False, False, False]},
                                   index=[5, 6, 8]),
-           ('x', 2): pd.DataFrame({'a': [7, 8, 9], 'b': [0, 0, 0]},
+           ('x', 2): pd.DataFrame({'a': [13094304034, 3489385935, 100006774],
+                                   'b': [0, 0, 0],
+                                   'c': [True, True, True]},
                                   index=[9, 9, 9])}
-    meta = make_meta({'a': 'i8', 'b': 'i8'}, index=pd.Index([], 'i8'))
+    meta = make_meta({'a': 'i8', 'b': 'i8', 'c': 'bool'}, index=pd.Index([], 'i8'))
     ddf1 = dd.DataFrame(dsk, 'x', meta, [0, 4, 9, 9])
     pdf1 = ddf1.compute()
 
@@ -627,9 +633,14 @@ def test_reductions(split_every):
     bools = pd.Series([True, False, True, False, True], dtype=bool)
     boolds = dd.from_pandas(bools, 2)
 
-    for dds, pds in [(ddf1.b, pdf1.b), (ddf1.a, pdf1.a),
-                     (ddf1['a'], pdf1['a']), (ddf1['b'], pdf1['b']),
-                     (nands1, nans1), (nands2, nans2), (nands3, nans3),
+    for dds, pds in [(ddf1.a, pdf1.a),
+                     (ddf1.b, pdf1.b),
+                     (ddf1.c, pdf1.c),
+                     (ddf1['a'], pdf1['a']),
+                     (ddf1['b'], pdf1['b']),
+                     (nands1, nans1),
+                     (nands2, nans2),
+                     (nands3, nans3),
                      (boolds, bools)]:
         assert isinstance(dds, dd.Series)
         assert isinstance(pds, pd.Series)
@@ -639,6 +650,7 @@ def test_reductions(split_every):
         assert_eq(dds.min(split_every=split_every), pds.min())
         assert_eq(dds.max(split_every=split_every), pds.max())
         assert_eq(dds.count(split_every=split_every), pds.count())
+
         with pytest.warns(None):
             # runtime warnings; https://github.com/dask/dask/issues/2381
             assert_eq(dds.std(split_every=split_every), pds.std())
@@ -648,6 +660,7 @@ def test_reductions(split_every):
         with pytest.warns(None):
             # runtime warnings; https://github.com/dask/dask/issues/2381
             assert_eq(dds.sem(split_every=split_every), pds.sem())
+
         assert_eq(dds.std(ddof=0, split_every=split_every), pds.std(ddof=0))
         assert_eq(dds.var(ddof=0, split_every=split_every), pds.var(ddof=0))
         assert_eq(dds.sem(ddof=0, split_every=split_every), pds.sem(ddof=0))
@@ -698,6 +711,23 @@ def test_reductions(split_every):
     assert_eq(ddf1.index.count(split_every=split_every), pd.notnull(pdf1.index).sum())
 
 
+@pytest.mark.parametrize('split_every', [False, 2])
+def test_reductions_timedelta(split_every):
+    ds = pd.Series(pd.to_timedelta([2, 3, 4, np.nan, 5]))
+    dds = dd.from_pandas(ds, 2)
+
+    with pytest.warns(None):
+        # runtime warnings; https://github.com/dask/dask/issues/2381
+        assert_eq(dds.var(split_every=split_every), ds.var())
+
+    assert_eq(dds.var(ddof=0, split_every=split_every), ds.var(ddof=0))
+
+    assert_eq(dds.sum(split_every=split_every), ds.sum())
+    assert_eq(dds.min(split_every=split_every), ds.min())
+    assert_eq(dds.max(split_every=split_every), ds.max())
+    assert_eq(dds.count(split_every=split_every), ds.count())
+
+
 @pytest.mark.parametrize('frame,axis,out',
                          [(pd.DataFrame({'a': [1, 2, 3],
                                          'b': [4, 5, 6]},
@@ -880,10 +910,9 @@ def check_raises(d, p, func):
     assert_eq(dds.max(), pds.max())
     assert_eq(dds.count(), pds.count())
 
-    # ToDo: pandas supports timedelta std, otherwise dask raises:
-    # incompatible type for a datetime/timedelta operation [__pow__]
+    # ToDo: pandas supports timedelta std, dask returns float64
     # assert_eq(dds.std(), pds.std())
-    # assert_eq(dds.var(), pds.var())
+    assert_eq(dds.var(), pds.var())
 
     # ToDo: pandas supports timedelta std, otherwise dask raises:
     # TypeError: unsupported operand type(s) for *: 'float' and 'Timedelta'
@@ -951,13 +980,22 @@ def test_reductions_frame(split_every):
     assert_dask_graph(ddf1.min(split_every=split_every), 'dataframe-min')
     assert_dask_graph(ddf1.max(split_every=split_every), 'dataframe-max')
     assert_dask_graph(ddf1.count(split_every=split_every), 'dataframe-count')
-    # std, var, sem, and mean consist of sum and count operations
-    assert_dask_graph(ddf1.std(split_every=split_every), 'dataframe-sum')
-    assert_dask_graph(ddf1.std(split_every=split_every), 'dataframe-count')
-    assert_dask_graph(ddf1.var(split_every=split_every), 'dataframe-sum')
-    assert_dask_graph(ddf1.var(split_every=split_every), 'dataframe-count')
-    assert_dask_graph(ddf1.sem(split_every=split_every), 'dataframe-sum')
-    assert_dask_graph(ddf1.sem(split_every=split_every), 'dataframe-count')
+
+    # std, var, sem, and mean consist of moment_* operations
+    assert_dask_graph(ddf1.std(split_every=split_every), 'dataframe-var')
+    assert_dask_graph(ddf1.std(split_every=split_every), 'moment_chunk')
+    assert_dask_graph(ddf1.std(split_every=split_every), 'moment_agg')
+    assert_dask_graph(ddf1.std(split_every=split_every), 'values')
+
+    assert_dask_graph(ddf1.var(split_every=split_every), 'moment_chunk')
+    assert_dask_graph(ddf1.var(split_every=split_every), 'moment_agg')
+    assert_dask_graph(ddf1.var(split_every=split_every), 'values')
+
+    assert_dask_graph(ddf1.sem(split_every=split_every), 'dataframe-var')
+    assert_dask_graph(ddf1.sem(split_every=split_every), 'moment_chunk')
+    assert_dask_graph(ddf1.sem(split_every=split_every), 'moment_agg')
+    assert_dask_graph(ddf1.sem(split_every=split_every), 'values')
+
     assert_dask_graph(ddf1.mean(split_every=split_every), 'dataframe-sum')
     assert_dask_graph(ddf1.mean(split_every=split_every), 'dataframe-count')
 
@@ -986,27 +1024,33 @@ def test_reductions_frame_dtypes():
     df = pd.DataFrame({'int': [1, 2, 3, 4, 5, 6, 7, 8],
                        'float': [1., 2., 3., 4., np.nan, 6., 7., 8.],
                        'dt': [pd.NaT] + [datetime(2011, i, 1) for i in range(1, 8)],
-                       'str': list('abcdefgh')})
+                       'str': list('abcdefgh'),
+                       'timedelta': pd.to_timedelta([1, 2, 3, 4, 5, 6, 7, 8]),
+                       'bool': [True, False] * 4
+                       })
 
     if HAS_INT_NA:
         df['intna'] = pd.array([1, 2, 3, 4, None, 6, 7, 8], dtype=pd.Int64Dtype())
 
     ddf = dd.from_pandas(df, 3)
+
+    # TODO: std and mean do not support timedelta dtype
+    df_no_timedelta = df.drop('timedelta', axis=1, inplace=False)
+    ddf_no_timedelta = dd.from_pandas(df_no_timedelta, 3)
+
     assert_eq(df.sum(), ddf.sum())
     assert_eq(df.prod(), ddf.prod())
     assert_eq(df.min(), ddf.min())
     assert_eq(df.max(), ddf.max())
     assert_eq(df.count(), ddf.count())
-    assert_eq(df.std(), ddf.std())
+    assert_eq(df_no_timedelta.std(), ddf_no_timedelta.std())
     assert_eq(df.var(), ddf.var())
     assert_eq(df.sem(), ddf.sem())
-    assert_eq(df.std(ddof=0), ddf.std(ddof=0))
+    assert_eq(df_no_timedelta.std(ddof=0), ddf_no_timedelta.std(ddof=0))
     assert_eq(df.var(ddof=0), ddf.var(ddof=0))
     assert_eq(df.sem(ddof=0), ddf.sem(ddof=0))
 
-    result = ddf.mean()
-    expected = df.mean()
-    assert_eq(expected, result)
+    assert_eq(df_no_timedelta.mean(), ddf_no_timedelta.mean())
 
     assert_eq(df._get_numeric_data(), ddf._get_numeric_data())
 
diff --git a/dask/dataframe/tests/test_dataframe.py b/dask/dataframe/tests/test_dataframe.py
index 0e18f5101ca..9beebdd0f1e 100644
--- a/dask/dataframe/tests/test_dataframe.py
+++ b/dask/dataframe/tests/test_dataframe.py
@@ -351,12 +351,6 @@ def test_describe(include, exclude, percentiles, subset):
     desc_ddf = ddf.describe(include=include, exclude=exclude, percentiles=percentiles)
     desc_df = df.describe(include=include, exclude=exclude, percentiles=percentiles)
 
-    # TODO: for timedelta columns there's overflow in var function, see #4233
-    # causing std to be nan, workaround this by dropping timedelta column before assertion
-    if 'f' in desc_ddf._meta:
-        desc_df = desc_df.drop('f', axis=1)
-        desc_ddf = desc_ddf.drop('f', axis=1)
-
     # Assert
     assert_eq(desc_ddf, desc_df)
 


Code review comments:
diff_hunk: @@ -1461,17 +1461,76 @@ def var(self, axis=None, skipna=True, ddof=1, split_every=False, dtype=None, out
                                     axis=axis, skipna=skipna, ddof=ddof)
             return handle_out(out, result)
         else:
-            num = self._get_numeric_data()
-            x = 1.0 * num.sum(skipna=skipna, split_every=split_every)
-            x2 = 1.0 * (num ** 2).sum(skipna=skipna, split_every=split_every)
-            n = num.count(split_every=split_every)
-            name = self._token_prefix + 'var'
-            result = map_partitions(methods.var_aggregate, x2, x, n,
-                                    token=name, meta=meta, ddof=ddof)
+            if self.ndim == 1:
+                result = self._var_1d(self, skipna, ddof, split_every)
+            elif len(self._meta_nonempty.select_dtypes(include=[np.timedelta64]).columns) > 0:
+                result = self._var_mixed(skipna, ddof, split_every)
+            else:
+                result = self._var_numeric(skipna, ddof, split_every)
+
             if isinstance(self, DataFrame):
                 result.divisions = (min(self.columns), max(self.columns))
             return handle_out(out, result)
 
+    def _var_numeric(self, skipna=True, ddof=1, split_every=False):
+        num = self.select_dtypes(include=['number', 'bool'])
+
+        values_dtype = num.values.dtype
+        array_values = num.values
+
+        if not np.issubdtype(values_dtype, np.number):
+            array_values = num.values.astype('f8')
+
+        var = da.nanvar if skipna or skipna is None else da.var
+        array_var = var(array_values, axis=0, ddof=ddof, split_every=split_every)
+
+        name = self._token_prefix + 'var-numeric' + tokenize(num, split_every)
+        cols = num._meta.columns if is_dataframe_like(num) else None
+
+        var_shape = num._meta_nonempty.values.var(axis=0).shape
+        array_var_name = (array_var._name,) + (0,) * len(var_shape)
+
+        layer = {(name, 0): (methods.wrap_var_reduction, array_var_name, cols)}
+        graph = HighLevelGraph.from_collections(name, layer, dependencies=[array_var])
+
+        return new_dd_object(graph, name, num._meta_nonempty.var(), divisions=[None, None])
+
+    def _var_mixed(self, skipna=True, ddof=1, split_every=False):
+        data = self.select_dtypes(include=['number', np.timedelta64, 'bool'])
+
+        vars = [self._var_1d(data[col_idx], skipna, ddof, split_every) for col_idx in data._meta.columns]
+        var_names = [(v._name, 0) for v in vars]
+
+        name = self._token_prefix + 'var-mixed-' + tokenize(data, split_every)
+
+        layer = {(name, 0): (methods.wrap_var_reduction, var_names, data._meta.columns)}
+        graph = HighLevelGraph.from_collections(name, layer, dependencies=vars)
+
+        return new_dd_object(graph, name, data._meta_nonempty.var(), divisions=[None, None])
+
+    def _var_1d(self, column, skipna=True, ddof=1, split_every=False):
+        is_timedelta = is_timedelta64_dtype(column._meta)
+
+        if is_timedelta:
+            column = column.dropna().astype('i8')
path: dask/dataframe/core.py
text: @user1:
Is this working around https://github.com/pandas-dev/pandas/issues/18880? Ideally we'd like to follow skipna, but that bug may prevent it.

@author:
@user1 looks like. I've looked into your PR for this bug and this what ideally should be done here as well:
```python
    if is_timedelta64_dtype(values) and not skipna:
        values = values.astype('float64')
        values[mask] = np.nan
```
But I suspect it's going to slow down things. What do you think?


@user1:
I just realized, we aren't using pandas' `nanstd`, so we don't need to worry about pandas-dev/pandas#18880 hopefully.

I'm not too concerned about the cost of computing the mask. Pandas does that internally as well, so it won't be outrageously slower (and this branch is only for `not skipna`, so it's not the default). Let's do something like that.

@author:
Done.