Fix machine error negative values on m2

[viraltux]

I found an issue with this test when m2 takes negative values due to machine error precision, see below:

x = [0.09999999999999981, 0.09999999999999978, 0.09999999999999978, 0.09999999999999978, 0.09999999999999977, 0.09999999999999983, 0.09999999999999981, 0.0999999999999998, 0.09999999999999981, 0.09999999999999978];

JarqueBeraTest(x)
ERROR: DomainError with -1.734723475976807e-18:
Exponentiation yielding a complex result requires a complex argument.
Replace x^y with (x+0im)^y, Complex(x)^y, or similar.

[devmotion]

I wonder what the (dis)advantages of this approach are compared with

Suggested change

	m2 = abs(-m1r^2 + m2r)
	m2 = max(-m1r^2 + m2r, 0)

[viraltux]

I was just about to pull a request for that option! It is better. At least mathematically more pleasant. Unless abs is faster than max... but probably not.

And now I am thinking that in the example that fails the second moment is not exactly zero which means that abs in that scenario is more precise than max that would set the second moment to exactly zero.

Now I think I should revert back to abs again! Okay, I'll sleep on it.

[viraltux]

Turns out that neither abs or max solves a NaN issue in skew and kurt when m2 equal zero.

The new change should improve performance since it calculates moments recursively replacing powers per products. In fact this approach prevents entirely negative numbers in the example that previously failed.

To deal with zeros in m2 the calculation stops and returns a default for zero JarqueBeraTest result.

[nalimilan]

The new change should improve performance since it calculates moments recursively replacing powers per products. In fact this approach prevents entirely negative numbers in the example that previously failed.

Have you tried only replacing yi^4 with yi^2 * yi^2? Julia has special cases for powers up to 3 (literal_pow) defined in terms of x*x*..., and I would expect the compiler to be able to optimize code automatically.

[nalimilan]

Also, could you add a test that fails without this PR?

[devmotion]

Have you tried only replacing yi^4 with yi^2 * yi^2?

Maybe (yi^2)^2 would be more performant? Or would it be more difficult for the compiler to optimize?

[nalimilan]

According to @code_native, these give the same result, so any of them is fine, but maybe yours is simpler to reason about for the compiler.

BTW, apparently precision will be better for x^4 than for other approaches, see JuliaLang/julia#20889 (comment). Not sure whether that's a problem here.