Maintenance will be performed on git.ligo.org, chat.ligo.org, containers.ligo.org, and docs.ligo.org on Tuesday 26 May 2020 starting at approximately 10am CDT. It is expected to take around 30 minutes and will involve a short period of downtime, around 5 minutes, towards the end of the maintenance period. Please address any questions, comments, or concerns to uwm-help@cgca.uwm.edu.

Commit 75944c67 authored by Leo Pound Singer's avatar Leo Pound Singer

Add sample code for finding credible regions

* Add sample code to test whether a sky position is in the 90%
  credible region.
* Add sample code to find the area of the 90% credible region.
parent 0b6f575f
......@@ -18,6 +18,10 @@ Version 9 (unreleased)
* Add tutorial on working with multi-resolution sky maps.
* Add sample code to test whether a sky position is in the 90% credible region.
* Add sample code to find the area of the 90% credible region.
Version 8 (2019-05-22)
----------------------
......
astropy-healpix
pygcn
healpy
ligo.skymap
seaborn
sphinx
git+https://github.com/djungelorm/sphinx-tabs@a2507e20b3ccae8d76adfcc39b4ca2515773ea19#egg=sphinx-tabs # https://github.com/djungelorm/sphinx-tabs/pull/43
......@@ -243,6 +243,73 @@ Let's find which pixel contains the point RA=194.95, Dec=27.98.
>>> ipix
13361492
Test if a Sky Location is in the 90% Credible Region
----------------------------------------------------
You can easily test if a given sky position is in the 90% credible region.
Let's continue using the sky position from the previous example, for which we
have already determined the pixel index.
Use the following simple algorithm to construct a map that gives the *credible
level of each pixel*:
1. Sort the pixels by descending probability density.
2. Cumulatively sum the pixels' probability.
3. Return the pixels to their original order.
In Python, you can use this simple recipe:
>>> i = np.flipud(np.argsort(hpx))
>>> sorted_credible_levels = np.cumsum(hpx[i])
>>> credible_levels = np.empty_like(sorted_credible_levels)
>>> credible_levels[i] = sorted_credible_levels
>>> credible_levels
array([1., 1., 1., ..., 1., 1., 1.])
.. note::
Observe that the values in the resulting *credible level map* vary inversely
with probability density: the most probable pixel is assigned to the
credible level 0.0, and the least likely pixel is assigned the credible
level 1.0.
.. tip::
This recipe is implemented in the package
:doc:`ligo.skymap <ligo.skymap:index>` as the function
:func:`~ligo.skymap.postprocess.util.find_greedy_credible_levels`:
>>> from ligo.skymap.postprocess import find_greedy_credible_levels
>>> credible_levels = find_greedy_credible_levels(hpx)
>>> credible_levels
array([1., 1., 1., ..., 1., 1., 1.])
To check if the pixel that we identified in the previous section is within the
90% credible level, simply test if the value of the credible level map is less
than or equal to 0.9 at that pixel:
>>> credible_levels[ipix]
0.9999999999947833
>>> credible_levels[ipix] <= 0.9
False
The credible level map has a value greater than 0.9 at that sky location,
therefore the sky location is outside the 90% credible region.
Find the Area of the 90% Credible Region
----------------------------------------
Since we just found the credible level map, it's easy to compute the 90%
credible area by counting the number of pixels inside the 90% credible region
and multiplying by the area per pixel.
In the Python expression below, note that ``(credible_levels <= 0.9)``
evaluates to a binary array; when it is summed over, true values are treated as
1 and false values are treated as 0.
>>> np.sum(credible_levels <= 0.9) * hp.nside2pixarea(nside, degrees=True)
30.979279207076633
Most Probable Sky Location
--------------------------
......
Markdown is supported
0%
or
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment