Poisson statistics with background - implementing wstat about sherpa

[cdeil]

This is something we need in gamma-ray astronomy as well.

@registerrier has implemented this and it's available here:
https://github.com/gammapy/gammapy/blob/master/gammapy/hspec/wstat.py
I don't know if it was @registerrier or someone else that filed issue #41.

It would be great if this could be implemented in Sherpa and we can just use it in Gammapy.
Is this something you plan to implement?
Or alternatively: can you post short instructions here how it should be implemented and one of us will start a pull request in the coming days / weeks?

[anetasie]

Yes, there is a plan to implement wstat in Sherpa. Great to see an implementation in Gammapy.

[dtnguyen2]

I implemented the wstat as described from the xspec site, and I'm currently
writing the C-extension (see file sherpa/include/stat_extension.hh). There
used to be three C-extension functions (staterrfct, statfct and
statfct_noerr) defined; However, only one of the functions has the comment
//PyNone MUST be incremented before being returned!!
and indeed the following code:
Py_INCREF(Py_None);

The question is, are the two other functions (statfct_noerr and statfct)
wrong for not incrementing the reference

PyObject* staterrfct( PyObject* self, PyObject* args )
{
....
// Py_None MUST be incremented before being returned!!
Py_INCREF(Py_None);
return Py_BuildValue( (char*)"(dO)", val, Py_None );

}

PyObject* statfct_noerr( PyObject* self, PyObject* args )
{
...
return Py_BuildValue( (char*)"(dO)", val, Py_None );
}

PyObject* statfct( PyObject* self, PyObject* args )
{
....
return Py_BuildValue( (char*)"(dN)", val, dev.return_new_ref() );

}

On Fri, Jul 31, 2015 at 12:48 PM, Aneta Siemiginowska <
[email protected]> wrote:

Yes, there is a plan to implement wstat in Sherpa. Great to see an
implementation in Gammapy.

—
Reply to this email directly or view it on GitHub
#73 (comment).

[DougBurke]

@dtnguyen2, in the version of stat_extension.hh I am looking at I see

• staterrfct template ends with https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/stat_extension.hh#L50

return err.return_new_ref();

• statfct_noerr ends with https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/stat_extension.hh#L103

// Py_None MUST be incremented before being returned!!
Py_INCREF(Py_None);
return Py_BuildValue( (char*)"(dO)", val, Py_None );

• statfct ends with https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/stat_extension.hh#L162

return Py_BuildValue( (char*)"(dN)", val, dev.return_new_ref() );

The return_new_ref method is defined in https://github.com/sherpa/sherpa/blob/master/sherpa/include/sherpa/array.hh#L75 and expands to a PyArray_Return statement, which creates a PyObject* and handles the reference counting for you.

The documentation for Py_BuildValue is rather dense and terse - https://docs.python.org/2/c-api/arg.html#c.Py_BuildValue - but the format string used here (of the form "(xx)") means that it returns a two-element tuple, with xx being either dO or dN. These mean that the first element is a double (and we don't have to bother with reference counting here as the value gets copied over to the object returned to Python), and the second is a PyObject*. The difference between O and N is that the former handles reference counting whereas the latter does not. As the return_new_ref method increases the reference count, this makes sense for statfct. For statfct_noerr it seems to make less sense (as O increases the reference count), but I am not an expert here. The introductory documentation https://docs.python.org/2/extending/extending.html#reference-counts doesn't really help, since it doesn't really consider anyone returning a tuple with one element being None. You can see from earlier examples - e.g. https://docs.python.org/2/extending/extending.html#back-to-the-example - that when returning Py_NONE directly it must have its reference count increased (or use the Py_RETURN_NONE macro). I'm guessing that this is the inspiration for the comment on line 103.

My hypothesis is that line 103 and 104 can be removed (so, in other words, the inverse of your question), but I have not tested it.

However, I'm not actually sure that the statfct_noerr template is ever used. When run from the top-level directory, I get:

% find sherpa -type f -exec grep STATFCT_NOERR {} /dev/null \;
sherpa/include/sherpa/stat_extension.hh:#define STATFCT_NOERR(name) _STATFCTSPEC(name, statfct_noerr) 

You can see in sherpa/stats/src/_statfcts.cc that statfct and staterrfct are used, but not the noerr variant.

[DougBurke]

I tried re-building master with STATFCT_NOERR and the statfct_noerr template commented out and the build completed, and all tests ran.

[DougBurke]

If you want to check for possible reference leaks, one thing I do is to run a small python program which calls the function you are checking repeatedly, and then just run top in another terminal and check that VIRT/RES/SHR/whatever they are called don't keep on increasing.

Here's an example program, which suggests that there's no problem with the staterrfct template:

% cat check_refs.py 
import numpy as np
from sherpa.stats import _statfcts
x = np.arange(1, 1000, 1)
while True:
    _statfcts.calc_chi2datavar_errors(x)

[dtnguyen2]

I'm more or less done with wstat. I just need a couple of things from
Aneta and/or Doug then I will commit:

1. The current error message says:
   raise FitErr( 'no bkg file is supplied, use cstat instead' )
   Let me know the exact message you want to tell the user.
2. Currently the doc string for the new WStat is very limited

class WStat(Likelihood):
"""
References
----------

.. [1] The description of the W statistic (`wstat`) in

https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html

"""

Let me know what you want me to put in there, check out the file
.../sherpa/stats/init.py for what the other stats class have etc...

On Thu, Aug 6, 2015 at 1:40 PM, Aneta Siemiginowska <
[email protected]> wrote:

Assigned #73 #73 to @dtnguyen2
https://github.com/dtnguyen2.

—
Reply to this email directly or view it on GitHub
#73 (comment).

[DougBurke]

@dtnguyen2 - how about https://gist.github.com/DougBurke/7fb922e782d82db921f7 as a start for WStat docs. I'm happy with the error message for the first case.

[dtnguyen2]

It took a while but I did find the mail where I was told that "I'm happy
with the error message for the first case". I don't have any problems
changing the error message but remember I did asked and was told that the
error message was fine. Here is the mail that I sent:

I'm more or less done with wstat. I just need a couple of things from
Aneta and/or Doug then I will commit:

1. The current error message says:
   raise FitErr( 'no bkg file is supplied, use cstat instead' )
   Let me know the exact message you want to tell the user.
2. Currently the doc string for the new WStat is very limited

class WStat(Likelihood):
"""
References
----------

.. [1] The description of the W statistic (`wstat`) in

https://heasarc.gsfc.nasa.gov/xanadu/xspec/manual/XSappendixStatistics.html

"""

Let me know what you want me to put in there, check out the file
.../sherpa/stats/init.py for what the other stats class have etc...

On Fri, Aug 21, 2015 at 2:45 PM, Doug Burke [email protected]
wrote:

@dtnguyen2 https://github.com/dtnguyen2 - how about
https://gist.github.com/DougBurke/7fb922e782d82db921f7 as a start for
WStat docs. I'm happy with the error message for the first case.

—
Reply to this email directly or view it on GitHub
#73 (comment).

[DougBurke]

Close by #94