fix cumulative_distribution_function

[EliasL]

Fixed a missing +1 to cumulative_distribution_function.

#111

[keflavich]

Thanks for the fix.

Could you add the example from #111 as a test? Your correct example appears to me to be correct, i.e., the final value in the CDF should be 1, so we should have a test that that condition is met.

I have to think through the searchsorted(..., right) logic more - your explanation makes sense, but I haven't touched this code in a decade and want to take a minute to think it over. Having a test case or two in place would really help.

[EliasL]

I'll see what i can do about the tests!
I can also offer an explanation for the right, left logic. Say we have data=[10, 20, 20, 50], and we can say these are scores from people playing a game. Then

CDF = searchsorted(data, data, side="left") / len(data)

gives [0. 0.25 0.25 0.75]. At each index, this is the percentage of people that got a score strictly lower than the score at this index. No one got a score lower than 10, so the first entry is 0.

CDF = searchsorted(data, data, side="right") / len(data)

gives [0.25 0.75 0.75 1. ], now showing the percentage of people that got a score lower than or equal to the score at this index. Three out of four people got a score lower than or equal to 20, so the second and third entry is 3/4.

Using "right" matches the docstring of the function:

Returns
    -------
    X : array
        The sorted, unique values in the data.
    probabilities : array
        The portion of the data that is less than or equal to X.

[EliasL]

I was just about to look into the tests when i found this comment in testing/test_powerlaw.py:

There is a subtle bug in the Clauset/plfit code involving the calculation of
the cumulative distribution function. Specifically, it assumes that only
discrete distributions can have repeat values and therefore performs the
calculation incorrectly in the case of a continuous distribution with repeated
values as occurs in the quakes and surnames data sets. The alpha values used
here for those data sets can be confirmed with the plfit code by forcing it use
the corresponding xmin values. Forcing powerlaw to calculate the cumulative
distribution function as done in the plfit code produces the same xmin values
as the plfit code and the other data sets produce identical results with both
plfit and powerlaw so the alpha and xmin values were changed as above for the
quakes and surnames data sets.

I don't know what this is about, but considering that it mentions the cumulative distribution function, whoever made these fixes may have fixed the wrong thing.

I have added a test to the cumulative_distribution_function (and the pdf for good measure).
The old CDF fails here:

# CDF gives probability of less than or equal to entry value.
# The last entry is the max of the sample. All other entries
# are therefore smaller than or equal, and the CDF
# should be 1.
assert_allclose(CDF[-1], 1.0)

But using the fixed CDF it now fails three other tests.

[keflavich]

Thanks for the explanations. It's entirely possible there were errors similar to what you've found here in the original plfit code that we incorrectly replicated.

I'm going to start a separate PR to get CI running so we can see these tests...

[keflavich]

Could you rebase against master? That will get CI running and will show the test failures.

I'd believe that you've corrected a long-standing error in the code, in which case we need to update the tests, but I don't want to jump to that conclusion without walking through the math myself or at least getting one other to review. I don't have time for that at the moment.

Curiously, the failures I see are in the lognormal and exponential tests. I don't see immediately why those fits would change so much.

The third test failure is this:

__________________________________________________________________________________________________________ FirstTestCase.test_power_law __________________________________________________________________________________________________________

self = <test_powerlaw.FirstTestCase testMethod=test_power_law>

    def test_power_law(self):
        print("Testing power law fits")

        rtol = .1
        atol = 0.01

        for k in references.keys():
            print(k)
            assert_allclose(results[k]['alpha'], references[k]['alpha'],
                            rtol=rtol, atol=atol, err_msg=k)

>           assert_allclose(results[k]['xmin'], references[k]['xmin'],
                            rtol=rtol, atol=atol, err_msg=k)
E           AssertionError:
E           Not equal to tolerance rtol=0.1, atol=0.01
E           words
E           Mismatched elements: 1 / 1 (100%)
E           Max absolute difference among violations: 3.
E           Max relative difference among violations: 0.42857143
E            ACTUAL: array(10.)
E            DESIRED: array(7)

testing/test_powerlaw.py:174: AssertionError

which looks like it might be exactly the off-by-one error for discrete (integer) valued data.

[Jfeatherstone]

Sorry it's taken me a while to get to this, I've been meaning to review it for a while!

I think the fix to the CDF calculation indeed does fix an old bug. I'm not incredibly familiar with the plfit package, but below is how I've convinced myself that the changes give better fitting. The example is written for v2.0.1 of powerlaw.

For the words dataset, as you've shown above, the old CDF calculation gives xmin=7, while the new CDF calculation gives xmin=10. The loglikelihood values are identical for both the old and new CDF calculation methods, and indicate that xmin=10 is actually the better fit.

data = powerlaw.load_test_dataset('words')

fit = powerlaw.Fit(data, xmin=10, discrete=True, estimate_discrete=False, test_all_xmin=True)
print(f'xmin=10: {fit.exponential.loglikelihood}')

fit = powerlaw.Fit(data, xmin=7, discrete=True, estimate_discrete=False, test_all_xmin=True)
print(f'xmin=7:  {fit.exponential.loglikelihood}')

(old CDF calculation)
xmin=10: -8974.760371952481
xmin=7: -11753.81757591329

(new CDF calculation)
xmin=10: -8974.760371952481
xmin=7: -11753.81757591329

I'll note that I initially started by looking at the Kolmogorov-Smirnov distance, which gives the following results:

(old CDF calculation)
xmin=10: 0.01185966250638315
xmin=7:  0.008256664404139191

(new CDF calculation)
xmin=10: 0.020094706055172384
xmin=7:  0.028951265966114303

ie. the opposite conclusion as above. That being said, since the KSD is dependent on the actual CDF calculation, it isn't a reliable way to calculate whether the CDF calculation itself is correct or not. The likelihood is calculated from the PDF instead, making it more suitable for this comparison.

This argument above only addresses the failure of the power law fitting test, though I think this exemplifies why caching the values of the loglikelihood ratios for distributions that don't actually fit the dataset isn't a good way to benchmark accuracy or consistency. I can't see any particular reason why the exponential or lognormal distributions fail while the others pass, other than just random chance.

---

Assuming we're convinced about this being a fix, there are a few things to do before merging this:

1. This PR needs to be rebased against master since v2.0 changed so much. I think since the changes are relatively compartmentalized, it might be easier just to hard reset to master and copy the new tests and the cdf() changes over. (The CDF changes are actually already included but commented out in v2.0+).
2. v2.0 already addresses the bounds on power law exponents, so we should remove those changes.
3. I think it would be a good idea to add some extra comments summarizing the discussion on left vs. right sorting and/or linking to this thread in the cdf() function.
4. How should we address the failing tests? One option would be to just increase the tolerances or even comment out those tests until they can be more permanently addressed. I have to make some changes to testing tolerances in a PR soon (mentioned in Saving and loading for Fit objects #126), so we could leave it until then.

@EliasL if you're busy and want me to make these changes, please let me know. Thanks for all of your work! :)

[EliasL]

Hello,

Your comments inspired me to make a propper test of the CDF, and this revealed that there was an even bigger issue in the KS distance itself.

I found a source which describes the KS distance more explicitly than what i have seen before:

$$\begin{aligned} D^{+} &= \max_i\left(\frac{i}{n} - U_{(i)}\right) \\\ D^{-} &= \max_i\left(U_{(i)} - \frac{i-1}{n}\right) \\\ D &= \max\left(D^{+}, D^{-}\right) \end{aligned}$$

Here we can recognize what the old code did in D- (i-1/n) and my new code in D+ ((i)/n). Note that they use one-indexing on i. It turns out that the real KS statistic uses both the old version (that i claimed was wrong), and my new proposal!

In the current KS code, D+ and D- are just doing (max(abs(CDF_diff)) with extra steps. Here is what is actually supposed to happen: When we have an empirical CDF (CDF_E), like for example [0.2, 0.4, 0.6, 0.8, 1.0], should the CDF_E be 0.2 or 0.0 at x=0.1? I have tried to illustrate the two options in this plot as CDF_E^+ and CDF_E^-:

The KS statistic suggests that we should try both, and keep whichever gives the largest distance.

There is an analytical expression for the expected KS distance from a model given a certain data set size n. With the help of AI, I find it to be

which works out to be

(To be confirmed...)

Using a sample size of n=20, we find that the expected KS distance (averaged over many samples) should be E=0.186, and that is exactly what we see with this new implementation of the KS distance.
I made a test which generates this plot: https://github.com/EliasL/powerlaw/blob/testing_old_CDF/testing/test_CDF.py

I know this is a lot, and I have not quite digested everything myself, but if you could take some time to see if you agree with me, I then suggest that we close this branch, and i'll fork the current head and try to suggest a fix to the KS function itself. (But I am not good with git, and open to any suggestions you might prefer)