ENH: numerically stable rolling_skew and rolling_kurt

[jaimefrio]

Close #6929. Hat tip to @behzadnouri for the z-scores idea in #8270, and to @seth-p for the general structure of the code, separating removal and addition of observations.

[jaimefrio]

Some rationale behind some of the design choices:

• I have joined roll_skew and roll_kurt into a single roll_higher_moment, since they shared most of the code.
• Since both skewness and kurtosis require a division by the variance, there is detection of exactly zero variance built into the function, to exactly detect NaNs.
• The actual computations of the skewness and kurtosis from the moments may look weird and inefficient, but they are trading extra multiplications for less divisions and/or square roots. This is the fastest I could manage, and doesn't seem to affect numerical stability.
• Rather than explicitly coding the four cases of val and prev being NaN or not, I have used @seth-p's idea of handling the removal and addition of an observation separately. This is slightly less efficient for the general case where neither is NaN, but I felt there was already enough complication. And the performance of the functions is either improved or maintained from current master, more details below. But there maya be a 15% of extra performance in there.
• Both kurtosis and skewness have terrible numerical stability when both the fourth or third central moment and the second, are very small. To minimize the effect on computations on some pathological cases, I have implemented @behzadnouri's idea of using z-scores as part of _rolling_func and _rolling_moment: they now take center_data and norm_data flags that trigger subtracting the mean and dividing by the standard deviation.
• To avoid unnecessarily slowing down execution, rolling_skew only centers the data, while rolling_kurt centers and normalizes it, see timings below.
• The loss of stability seems to be related with the removal of observations, so neither expanding_skew nor expanding_kurt do any centering or normalizing, although it would be very easy to implement.

Timings below are on x = np.random.rand(1e6) calling rolling_xxx(x, 100, 0). The chosen implementation is marked with ** **. Performance is unchanged for rolling_skew and almost:doubled for rolling_kurt:

                rolling_skew rolling_kurt
master              55 ms       129 ms
data as-is          50 ms        47 ms
data centered     **56 ms**      54 ms
data normalized     71 ms      **69 ms**

[seth-p]

Would it make sense to an a ''ddof'' argument?

On Sep 29, 2014, at 11:48 PM, Jaime notifications@github.com wrote:

Some rationale behind some of the design choices:

I have joined roll_skew and roll_kurt into a single roll_higher_moment, since they shared most of the code.
Since both skewness and kurtosis require a division by the variance, there is detection of exactly zero variance built into the function, to exactly detect NaNs.
The actual computations of the skewness and kurtosis from the moments may look weird and inefficient, but they are trading extra multiplications for less divisions and/or square roots. This is the fastest I could manage, and doesn't seem to affect numerical stability.
Rather than explicitly coding the four cases of val and prev being NaN or not, I have used @seth-p's idea of handling the removal and addition of an observation separately. This is slightly less efficient for the general case where neither is NaN, but I felt there was already enough complication. And the performance of the functions is either improved or maintained from current master, more details below. But there maya be a 15% of extra performance in there.
Both kurtosis and skewness have terrible numerical stability when both the fourth or third central moment and the second, are very small. To minimize the effect on computations on some pathological cases, I have implemented @behzadnouri's idea of using z-scores as part of _rolling_func and _rolling_moment: they now take center_data and norm_data flags that trigger subtracting the mean and dividing by the standard deviation.
To avoid unnecessarily slowing down execution, rolling_skew only centers the data, while rolling_kurt centers and normalizes it, see timings below.
The loss of stability seems to be related with the removal of observations, so neither expanding_skew nor expanding_kurt do any centering or normalizing, although it would be very easy to implement.
Timings below are on x = np.random.rand(1e6) calling rolling_xxx(x, 100, 0). The chosen implementation is marked with ** **. Performance is unchanged for rolling_skew and almost:doubled for rolling_kurt:

            rolling_skew rolling_kurt

master 55 ms 129 ms
data as-is 50 ms 47 ms
data centered 56 ms 54 ms
data normalized 71 ms 69 ms
—
Reply to this email directly or view it on GitHub.

[seth-p]

Also, I'd suggest adding consistency checks on rolling/expanding_skew/kurt (both ddof=0 and ddof=1) to _test_moments_consistency(), test_expanding_consistency(), and test_rolling_consistency().

[jaimefrio]

I don't think it is well defined what ddof would mean here... We could add a bias=True/False following scipy, but everything in pandas is unbiased by default, not sure if it would make much sense.

[seth-p]

I'm not sure I agree that Inf values should be converted to NaN. I'm not sure they shouldn't, but my initial reaction is that they shouldn't be, ie should be left as is.

Also, rather than explicitly masking out NaN values, why not use np.nanmean() and np.nanstd()? (I don't recall in which version of numpy they were introduced.) I suspect that will be quite a bit faster.

[seth-p]

Everything is unbiased by default, but (now, as of 0.15.0) is specifyable. Seems like these should be too.

[jaimefrio]

np.nanmean and np.nanstd were included in numpy 1.8, but they mostly do the same as here, you can look at the source here. And if we called nanmean and nanstd we would be creating (almost) the same mask three times, which is what I was trying to avoid.

Also, if you don't mask out the np.infs they mess up the rolling-ness of the calculation, and everything will be np.inf or np.nan thereafter. This was already in there, I just refactored it to reuse the mask. It can be changed, but we should change it for all the rolling functions, not just for these two: I don;t think it belongs in this PR.

[jreback]

put the defaults here (after bool),e .g. center_data : bool, default False

[seth-p]

@jaimefrio, I'm not sure how ddof would figure into the calculation of skewness and kurtosis, but I think it would be easy to support a bias parameter analogous to that in ewmvar(). Following the notation in the second page of http://www.jstor.org/discover/10.2307/2988433 (page "184), I think the existing (bias=False) behavior is to calculate skew and kurtosis as G1 = K3 / K2^(3/2) and G2 = K4 / (K2^2), respectively. I'd then suggest that the bias=True behavior would be to calculate skew and kurtosis as m3 / m2^(3/2) and m4 / (m2^2) - 3, respectively. Though perhaps bias=True isn't the right way to refer to the latter, as m3 / m2^(3/2) isn't actually biased. At any rate, with these "biased" estimators you could then have simple consistency tests of moments on the lines of the existing check that biased var(x) == mean(x^2) - mean(x)^2 in _test_moments_consistency():

            mean_x = mean(x)
            ...
            for (std, var, cov, skew, kurt) in [
                (std_biased, var_biased, cov_biased, skew_biased, kurt_biased),
                (std_unbiased, var_unbiased, cov_unbiased, skew_unbiased, kurt_unbiased)]:
                var_x = var(x)
                std_x = std(x)
                skew_x = skew(x)
                kurt_x = kurt(x)

                ...
                if var is var_biased:
                    # check that biased var(x) == mean(x^2) - mean(x)^2
                    mean_x2 = mean(x * x)
                    assert_equal(var_x, mean_x2 - (mean_x * mean_x))
                    # check that "biased" skew(x) == third_moment / second_moment^(3/2)
                    third_moment_x = mean(x**3) - 3 * mean_x2 * mean_x + 3 * mean_x**3 - mean_x**3 # should collapse last two terms to + 2 * mean_x**3
                    assert_equal(skew_x, third_moment_x / (std_x**3))
                    # check that "biased" kurt(x) == fourth_moment / second_moment^2 - 3
                    fourth_moment_x = mean(x**4) - 4 * mean(x**3) * mean_x + 6 * mean_x2 * mean_x**2 - 4 * mean_x**4 + mean_x**4 # should collapse last two terms to - 3 * mean_x**4
                    assert_equal(skew_x, fourth_moment / var_x**2 - 3.)

There are probably typos in the above code, so caveat emptor.

[jaimefrio]

I'll just replicate what scipy.stats does, which is mostly what you describe.

[jreback]

@jaimefrio can you rebase/revisit ?

[jreback]

@jaimefrio can you rebase and add a release note, lgtm otherwise.

[jreback]

@jaimefrio can you rebase / revist ?

[jaimefrio]

Yes, will do, @jreback. Sorry I missed your earlier pings!

[jreback]

@jaimefrio can you update for 0.17.0

[jreback]

can you update?

[jreback]

@jaimefrio can you rebase?

[jreback]

@jaimefrio can you rebase

[jaimefrio]

I've rebased, but I don't remember anymore what was missing for this to be ready...

[jreback]

can you add a versionadded directiver here

[jreback]

can you add the examples from top of the issue #6929 e.g. validating the problem skew/kurt results that we were getting before (IOW, calculate and compare to what they should be)

[jreback]

can you add a release note in v0.17.0 (maybe a separate mini-section), to show what has changed? (e.g. you can use the same examples that show the stability problem).

[jreback]

@jaimefrio can you update

[jreback]

can you add a release note to 0.17.1

[jreback]

@jaimefrio can you rebase/update?

[jreback]

closing, but pls reopen when you can update