Circstats issue3705 - #4472

Thank for another pull request. Some of these functions are already in scipy and I do not think we want to duplicate the work here unless we make the functions specifically useful for astronomy/astropy. So I think the inclusion of these statistics does require some thoughtful discussion.

One option would be that you should submit this pull request to scipy, which would make available to the large community and probably a good idea, as scipy already as a circstats module.

However, a useful thing though would be making sure the stat functions work with astropy.qunaitity and that might be a good reason for their inclusion here. In addition, you have added some more functionality which is not available in the scipy.stats.

Also, are you considering implementing more of the circularly functionality from R?

Thanks @crawfordsm for commenting.

In fact, the circmean function is pretty much the same as the scipy one. However, the circvar function is not quite a duplication of the scipy one since they use a different approach to define "dispersion" for circular data.

I will definitely see what is necessary to the functions so that they work properly with astropy.quantity.

Yes. Do you have any circular statistics functionality in mind?

Just some nitpicking:

• Minor PEP8 formatting stuff.
• Perhaps it is better to replace AssertionError with ValueError and more descriptive error messages.
• The documentation is very nice but all buried within API docstring. They should be moved to a proper RST file as part of stats documentation.

@astronomeralex, @eteq, and @mhvk -- I'm curious about your thoughts about inclusion of the code:

pros: It is code useful for pulsar astronomy as mentioned in #3705 and has several functions that are not currently available elsewhere.

cons: some functionality overlaps with scipy.stats and it is more generally statistical methods rather than astronomy specific.

@crawfordsm @eteq @mhvk That's a good question. I would lean towards including the code for your pro's listed above. This, and other pull requests, keep pressing this issue of whether stats are astro-specific. It seems prudent that we solve this question for the whole module instead of pull by pull.

@crawfordsm - I agree with you: these routines are useful outside of astronomy as well and the aim should be to have them in scipy. The main objection I can see is that for some it would be useful to have them available now rather than wait for things to have moved to a stable release in scipy. For that, we could have our versions in a utils.compat.scipy module, in analogy with what we currently do for numpy (see #2327).

In the process, we would hopefully also resolve any issues with scipy not being compatible with Quantity (and therefore likely having problems with other ndarray subclasses), which would benefit everyone as well. Indeed, @mirca, if you adapt your tests to work with Quantity, I'd love to have those as part of tests for astropy.units!

@mirca - is this something you are willing to pursue? At some level, you'd be our guinea pig, to see how such PRs fare.

p.s. For the same reasons as described above, I do think it would be good to eventually move astropy.units to scipy as well. My main hesitation so far has been that, without __numpy_ufunc__, the state of Quantity is not yet as ideal as it could be.

@mirca - yes, one would like to be sure the code works with Quantity, including if the quantity is in degrees. Since the first thing the present code does is turn any input into a straight ndarray, this is unlikely to work as is. Note that I think that here one should just ducktype (or at least do np.asanyarray); most of the arithmatic should just work with quantities.

OK, I did what I should have done before commenting too much. First, I tried and the scipy routine definitely is not Quantity safe. Second, I looked at the scipy code, and think it will be somewhat problematic to make this work with Quantity, mostly because they use high and low values that have defaults that are float rather than of some unit, and are clearly meant to do the unit conversion that is automatic if one uses Quantity. I worry though that they would be less than keen to start inserting that are quantity-specific.

Given the above and given that the code they use is relatively trivial and even easier to write using Quantity, it may just makes sense to do that after all.

That said, I do think any code in astropy should fully use Quantity. I'm also wary of using definitions that are different from those used in scipy (even if they are consistent with what is used in R).

Just to show how the code might look like if full Quantity-aware

def circmean(sample, axis=None):
    sample = Longitude(sample)
    with u.add_enabled_equivalencies(u.dimensionless_angles()):
        complex_mean = np.exp(1j * sample).mean(axis=axis)
    return Longitude(np.arctan2(complex_mean.imag, complex_mean.real),
                     sample.unit, wrap_angle=sample.wrap_angle)

from astropy.coordinates import Longitude
theta = Longitude(np.random.uniform(np.pi,-np.pi,10), u.rad)
circmean(theta)
# <Longitude 4.086819352329046 rad>
circmean(theta.to(u.deg))
# <Longitude 234.15750052084297 deg>
theta.to(u.deg)
# <Longitude [ 318.40807173, 258.09501963, 205.94616268, 256.10717963,
#               86.79577277,  19.49205343, 186.6603956 , 231.79219607,
#              109.07587581, 249.56135721] deg>

Alternatively, one could wrap the scipy code where possible:

def circmean(sample, axis=None):
    sample = Longitude(sample)
    circmean = scipy.stats.circmean(sample.to(u.rad.value), axis=axis)
    return Longitude(circmean * u.rad, sample.unit, wrap_angle=sample.wrap_angle)

@mirca Once you have updated the code to work with Quantities, I'll give it a final review before merging. Although it looks like you will need to rebase most likely due to updates to the change log.

Hi @crawfordsm, @mhvk. I've made some modifications and now the functions are able to handle both Quantity and non-Quantity types. More precisely, if the input data is a Quantity object, then the returned value, when appropriate, will be also a Quantity object. On the other hand, if the input data is not a Quantity object, then the returned value will not be a Quantity either. Is that the desired behavior?
I also add a parameter called deg, in which the user should set to True whenever the circular data is expressed in degrees.

@mirca - had only the most cursory of looks so far, and liked very much that for most of the code there is, as I would have hoped, absolutely no difference between what is done for ndarray and for Quantity. This is duck-typing at its best!

The exceptions are the degree to radian conversions; my sense would be to not allow this. Effectively, this is equivalent to allowing some special units to be attached to ndarray, which is just what Quantity is for. So, while I'm quite happy for plain numbers in an ndarray to be considered to have dimension of radian (as this is equivalent to dimensionless, which is what a plain number should be), I don't like the idea of allowing plain numbers to really be degrees.

@crawfordsm, @mhvk, Does this look better? (:

I made some modifications in such a way that all functions now will return Quantity regardless whether the user's input is a ndarray or a Quantity object. That basically removes the need to convert deg to rad.

@mirca I will try to review it this week and give you comments.

@crawfordsm any thoughts?

Thanks for all the updates @mirca.

A final note is that the docstrings need to be made consistent with how the functions work. I've left an example for circmean, but could you update all functions where it is appropriate? Other than that, I'll be happy to have this merged.

@crawfordsm Nice catch! Do I need to rebase?

Looks like it

@crawfordsm done! :)

Thanks @mirca !