QU-fitting PA uncertainty does not take angle wrapping into account
jackieykma opened this issue · 5 comments
While the current implementation of QU-fitting takes into account of PA wrapping during the sampling (with, e.g., nested sampling), I suspect that it does not take the angle wrapping into account while calculating the uncertainties of PA
I attach an example source Stokes qu (i.e. fractional) spectrum here as J182900-002018.txt. Upon running QU-fitting using a Faraday-thin model (m1), the resulting PA is 11.26 +164.71 - 7.81 deg. Looking at the corner plot here J182900-002018_ciradafig_m1_corner.pdf, one would expect the uncertainty to be rather symmetric on both sides, if the angle wrapping has been taken into account
I'm not sure yet whether this is an issue within the RM-Tools, or within the sampler. I used Nestle (v 0.2.0) through bilby (v 1.2.1)
Hi @jackieykma - this is interesting. I think we could probably take this into account within RM-Tools. EDIT: but it would be nice for the sampler to just handle it.
I did just find this in the bilby docs:
https://lscsoft.docs.ligo.org/bilby/api/bilby.core.prior.base.Constraint.html
The boundary condition of the prior, can be ‘periodic’, ‘reflective’ Currently implemented in cpnest, dynesty and pymultinest.
@jackieykma would you mind re-running this example with the pymultinest and/or dynesty samplers?
Thanks @AlecThomson ! See results below (also including nestle for reference)
J182900-002018fig_m1_corner.pdf
J182900-002018fig_m1_corner.pdf
J182900-002018fig_m1_corner.pdf
What I notice is that:
- All have them seem to experience the issue of vastly overestimated PA error
- (Kind of off-topic) The results from nestle & dynesty agree with each other well, while those from pymultinest deviate slightly from the other two
Following up on this issue, I've looked through the QU-fitting code and found that the reported values & uncertainties for each parameter is calculated in the following line
RM-Tools/RMtools_1D/do_QUfit_1D_mnest.py
Line 429 in a52419a
This means that the calculations (which does not take angle wrapping into account) stem from the get_one_dimensional_median_and_error_bar method of bilby
I am not sure if it is worth contacting the bilby people about this issue
Instead, I propose that we can slightly modify the QU-fitting module of RM-Tools to get around this
Specifically, we can obtain two sets of uncertainty values instead of one
- One directly from the sample, as we do now
- Another one from shifting the sample PA values by 90 deg
Out of these two sets of uncertainties, we can take the minimum out of the two
This would solve the issue because, for sources with PA near the 0 / 180 deg boundary, the second sample will be shifted to around 90 deg (where there would be no angle wrapping)
I'm happy to dig into this and create a pull request
I think the idea of a 90° rotation and re-calculation is very good. I think there's a few questions that would have to be sorted out:
- at what stage can we do this (in order to minimize re-calculation)? Definitely not re-doing the whole nested sampling. But if it can be done to the samples, that should be pretty efficient?
- Is there an easy way to identify which parameters need this to be done? Some models have more than one angle parameter, so this solution requires finding the necessary parameters to do this to.
Both problems are pretty easy, I suspect, but I don't know this part of the code very well so I can't say exactly what changes would be required.
If you're happy to do it, I'm happy to let you do it. 😁
Hi @Cameron-Van-Eck ,
- Agreed --- I don't think we should re-do the nested sampling, especially since we have found that QU-fitting can take a long time to run for the more complex models. Editing the samples is probably the way to go
- Yes --- I try to find parameters with name ending with "deg" (which is the units of the parameter), and only do the 90 deg shift for those
Just now I have coded it up and seems to work fine on one source from my VLA data (see attached for before v.s. after)
I will test this further with simulated sources before submitting a pull request