nanograv / enterprise

ENTERPRISE (Enhanced Numerical Toolbox Enabling a Robust PulsaR Inference SuitE) is a pulsar timing analysis code, aimed at noise analysis, gravitational-wave searches, and timing model analysis.
https://enterprise.readthedocs.io
MIT License
64 stars 65 forks source link

How to do wideband TOAs #129

Open jellis18 opened 6 years ago

jellis18 commented 6 years ago

I've started to think of how we are going to do wideband TOAs. In this case we have both TOAs and DM measurements per band per epoch. In PAL2 I just appended the extra DM measurements onto the residuals and constructed the components as follows

eq

where d is the data vector (dt and dx are the time residuals and DM residuals, respectively), T is the overall basis matrix with T' the standard TOA basis matrix (not including DMX), K is the DMX design matrix and K is the DMX design matrix again but without the 1/nu^2 scaling and N is the new white noise covariance matrix which is just the normal TOA uncertainties with the DM uncertainties appended.

Once I had this, the likelihood is the same as before since everything is just a product of these three components. Since the signal is built up from multiple components instead of defining T "by hand" as was done in PAL2 this approach is a bit harder.

I'm not quite sure of an elegant way to do this. Any suggestions?

vallis commented 6 years ago

Questions....

jellis18 commented 6 years ago

In the above you define K twice... I think you meant to say that D is the DMX design matrix without frequency scaling... but shouldn't that be just the identity matrix?

Yes I meant to say D but it isn't the identity matrix but the an exploder type matrix. The wideband TOAs still have one point per band per epoch. So if we have L-band and 800 MHz data then each epoch will have two TOAs and two DMs.

Are TOAs and DMs really uncorrelated? Is that achieved by referring each DM to a different frequency, and if so, are we rescaling appropriately in DMX?

I've never quite understood this but yes, apparently they are constructed so that they are uncorrelated. I don't think you would need to do anything to DMX..

If the DM estimate is already reflected in the TOAs, can't we handle the DMs (and marginalize over them) simply as timing-model parameters with Gaussian priors? There would be one such parameter per epoch, but that's OK...

This is what we thought of doing originally but the problem becomes that per epoch we only have 1 DMX parameter but we have two DM measurements (assuming 2 TOAs per epoch) so if the DMX model in the TOAs is Ku where K is the DMX basis and u are the actual DMX parameters, normally we would want a prior on u (i.e. a mean and variance per epoch) but instead we have two DM measurements per epoch (one for each band)

vallis commented 6 years ago

If the two DM measurements are truly uncorrelated with the TOAs and if we believe their normal errors, then we would be able to take their error-weighted average as the central value of the DMX parameter, and their combined error as the width of the prior. (Our fancier math should achieve the same)

In practice, we need to be able to handle timing-model parameters with a mean value other than 0, and we need to understand what DM effect, if any is subtracted from the TOA that we get in the dataset.

What's a good reference for the wideband timing?

vallis commented 6 years ago

Another way of thinking about it is that each DMX parameter gets a combined prior given by all of its measurements at that epoch.

jellis18 commented 6 years ago

Here is the reference.

demorest commented 6 years ago

Hi guys, the uncorrelated TOA and DM is explained in Tim's paper (http://adsabs.harvard.edu/abs/2014ApJ...790...93P). It's just what Michele said, you can choose a TOA reference frequency that results in no correlation between TOA and DM.

I think the DM uncertainties are probably OK but there are still some issues we don't fully understand about systematic biases between DMs measured in the different frequency bands. We might want to be able to include something like "DMJUMP" in the model to account for this. Not sure if this would work in the "average DMs to get a DMX prior" approach?

vallis commented 6 years ago

Thanks Paul. That means that every TOA now carries the reference frequency (under which .tim flag?), and we must use that to apply further DM corrections.

What would a DMJUMP parameter do exactly?

jellis18 commented 6 years ago

That means that every TOA now carries the reference frequency (under which .tim flag?), and we must use that to apply further DM corrections.

As far as I know, nothing is done with the reference frequencies unless internally somehow

What would a DMJUMP parameter do exactly?

Maybe this isn't what you were asking but a DMJUMP basically just fits a mean per frequency band, thus removing systematic offsets.

Overall, to me, the way I mentioned doing this in the first comment of the thread is the "correct" way to do this (maybe with other components added like DMJUMP, I already use a DMEFAC) in that we basically have a vector valued input now, a TOA and DM per "observation" and the model for the DM is

dx = Du + n_dm

and the model for the TOAs is

dt = T'b + Ku + n_toa

thus the DM is completely described by the "DMX" model but the TOAs have other delays as well. I could be wrong about this and there may be a better way of doing this...

The problem is I don't know a nice way of getting this in the enterprise framework.

demorest commented 6 years ago

Every TOA has a frequency associated with it. It's not a flag, it's one of the required pieces of info (it's the 2nd field in tempo2-format TOA lines). For standard TOAs this is usually set to the center freq of whatever band was averaged to make the TOA. For wide-band TOAs we set it equal to the DM reference freq.

Also, it's been a little while since I thought this through, but (if we temporarily ignore complications like DMJUMP, DMEFAC, etc) doesn't using the measured DMs as priors on DMX turn out to be mathematically equivalent to the "TOA+DM fit" approach Justin wrote up above?

jellis18 commented 6 years ago

@demorest, I'm sorry that I can never seem to understand this but for DMX there is one value per epoch and for the measured DMs there are two values per epoch. How can you use the two values as a prior? Is it

p(DMX|DM) ~ \prod_{i=0}^{i=2} \exp(-1/2 (DMX-DM_i)^2/sigma_DM_i^2)

demorest commented 6 years ago

I think it would just be what Michele said, take the appropriate weighted average of all measured DM values available within the DMX bin.

vallis commented 6 years ago

Yes! After all, one man's likelihood (of DMX given the DM measurements) is another man's likelihood. The key point is that the DM measurement do not depend on the TOAs. Or at least that's what is claimed!

In terms of implementation in Enterprise, though, while we can populate the combined sigma in the Phi matrix for the timing model, we do need to subtract the combined mean, multiplied by the design matrix column (really two elements). Perhaps that can be handled by the timing-model class.

Now, if you add a MEANDM_i as a timing-model parameter for each band i, the accounting gets more complicated—specifically, the prior covariance matrix for the timing-model parameters becomes non-diagonal.