ENH: stats.cauchy: improve cdf and quantile accuracy in tail about scipy

[WarrenWeckesser]

I found that the standard form (i.e. loc=0, scale=1) of the CDF and SF for the Cauchy distribution can be written as follows:

def cauchy_cdf(x):
    return np.arctan2(1, -x)/np.pi


def cauchy_sf(x):
    return np.arctan2(1, x)/np.pi

These give good precision over the entire domain.

The inverses can be implemented by wrapping the functions available in Boost.

I'll submit a PR soon with these updates.

[tk-yoshimura]

For example, how about the following code:

    def _cdf(self, x): 
         return np.where(
             x > -2.0, 
             0.5 + np.arctan(x) / np.pi, 
             -np.arctan(1.0 / x) / np.pi
         )
    
    def _sf(self, x): 
        return np.where(
             x < 2.0, 
             0.5 - np.arctan(x) / np.pi, 
             np.arctan(1.0 / x) / np.pi
        )
    
    def _ppf(self, q): 
        return -1.0 / np.tan(q * np.pi)
    
    def _isf(self, q): 
        return 1.0 / np.tan(q * np.pi)

I wish numpy had integer periodic trigonometric functions such as sinpi as implemented in MATLAB.

[dschmitz89]

Thanks @tk-yoshimura , this looks like a good suggestion. One minor thing is to use scipy's lazywhere instead of np.where (see for example here) for better performance. Can we interest you in a PR? We welcome expertise on precision related issues.

[tk-yoshimura]

Thanks @dschmitz89, you have also been helpful with the yule-simon distribution, non-central t distribution, and zipf distribution.
I would also like to help scipy with the stable distribution I implemented in C#, so I will start by forking to understand the structure of scipy.

It may take some time to issue the PR, so if you are in a hurry, may I ask for your help?
I am better at C++,C# than python.

If _lazywhere is used, it will look like this:
Is this correct?

    def _cdf(self, x): 
         return _lazywhere(
             x > -2.0, [x],
             lambda x: 0.5 + np.arctan(x) / np.pi, 
             f2=lambda x: -np.arctan(1.0 / x) / np.pi
         )
    
    def _sf(self, x): 
        return _lazywhere(
             x < 2.0, [x],
             lambda x: 0.5 - np.arctan(x) / np.pi, 
             f2=lambda x: np.arctan(1.0 / x) / np.pi
        )

[dschmitz89]

Thanks @dschmitz89, you have also been helpful with the yule-simon distribution, non-central t distribution, and zipf distribution. I would also like to help scipy with the stable distribution I implemented in C#, so I will start by forking to understand the structure of scipy.

It may take some time to issue the PR, so if you are in a hurry, may I ask for your help? I am better at C++,C# than python.

If _lazywhere is used, it will look like this: Is this correct?

    def _cdf(self, x): 
         return _lazywhere(
             x > -2.0, [x],
             lambda x: 0.5 + np.arctan(x) / np.pi, 
             f2=lambda x: -np.arctan(1.0 / x) / np.pi
         )
    
    def _sf(self, x): 
        return _lazywhere(
             x < 2.0, [x],
             lambda x: 0.5 - np.arctan(x) / np.pi, 
             f2=lambda x: np.arctan(1.0 / x) / np.pi
        )

We are not in a hurry, don't worry. Besides quick fixes we also seek to expand our contributor base :).

I think the code snippet should look like:

    def _cdf(self, x): 
         return _lazywhere(
             x > -2.0, (x, ),
             lambda x: 0.5 + np.arctan(x) / np.pi, 
             f2=lambda x: -np.arctan(1.0 / x) / np.pi
         )
    
    def _sf(self, x): 
        return _lazywhere(
             x < 2.0, (x, ),
             lambda x: 0.5 - np.arctan(x) / np.pi, 
             f2=lambda x: np.arctan(1.0 / x) / np.pi
        )

So you were almost correct, only (x, ) instead of [x] is required.

For your PR please also include a test that currently fails and passes with your update in this file.

[nickodell]

@dschmitz89

Just for my own curiosity, why is it better to use a tuple vs. a list here? Is it a speed thing? Is it better style?

[fancidev]

I wish numpy had integer periodic trigonometric functions such as sinpi as implemented in MATLAB.

SciPy does have sindg and related functions, taking degrees as arguments. Though under the hood they are converted to radian (after range reduction) and fed into the usual trigonometric functions.

For the functions here range reduction is not an issue since 0 <= q <= 1. Apart from producing certain special values exactly, such as tanPi(0.25)==1, are there other benefits of say tanPi over tan?

[tk-yoshimura]

@fancidev

The problem with my code is that the following is not exactly zero.
The reason for this is that 0.5*np.pi is not exactly pi/2 in floating point numbers.

ppf(0.5)=-6.123233995736766e-17

[fancidev]

The problem with my code is that the following is not exactly zero. The reason for this is that 0.5*np.pi is not exactly pi/2 in floating point numbers.

ppf(0.5)=-6.123233995736766e-17

I see. I feel it would not be a bad idea to add tanPi and related functions to SciPy, but I’m not sure if it would be attractive to the community. Maybe I should ask that in the mailing list.

Update: I just did. https://discuss.scientific-python.org/t/rfc-add-sinpi-cospi-tanpi-and-cotpi-functions-to-special/1193

[dschmitz89]

@dschmitz89

Just for my own curiosity, why is it better to use a tuple vs. a list here? Is it a speed thing? Is it better style?

@nickodell I think at some point I had tried it in the past and it did not work. I am not 100% sure anymore though, so feel free to test a list instead of the tuple.