Chromatic Gaussian event model

[jeremy-baier]

Adds a transient chromatic Gaussian event model as a delay signal. See https://iopscience.iop.org/article/10.3847/2041-8213/acdd03/pdf for details on the model, but note a typo missing the negative sign in the definition of the Gaussian.

Derivatives checked here:
https://www.wolframcloud.com/obj/f6b351f0-5c17-4777-ab4e-fd0f6ff03c41

[jeremy-baier]

@davidkaplantest
I think I have a working draft of this now. Would you be able to review this?
I have written tests and gotten those to pass.
I would appreciate it if you could check my math on the log10sigma derivative. I got a bit confused there.

[dlakaplan]

I'll start taking a peak at this. Note that the docs render well:
https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent
but I think you could probably use more standard formatting for the paper reference (at least a URL).

[dlakaplan]

I see differences in the parameters and the analytical form (no sign in the form, amplitude is not log10). Can you connect the analytic expression with the parameters more explicitly?

[dlakaplan]

@davidkaplantest I think I have a working draft of this now. Would you be able to review this? I have written tests and gotten those to pass. I would appreciate it if you could check my math on the log10sigma derivative. I got a bit confused there.

I can look at this, but I'd say that if you've gotten the fitting to work then that's a good sign that the derivatives are correct (or at least close enough).

[jeremy-baier]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.

I will plan to rewrite the docs so that it matched the parameters themselves.

[dlakaplan]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.

I will plan to rewrite the docs so that it matched the parameters themselves.

Great. Note that it can be frustrating to check the builds since it takes an ~hour for GitHub to do it. If you really want to dig into that I can give you instructions for doing a local build which can help iterate faster.

[jeremy-baier]

I finally sat down and ran a simple test of this. It seems that the model recovery happens very well with a toy model.
The attached images show 1) injected chromatic Gaussian 2) intentionally slightly misspecified timing model with injected chromatic Gaussian 3).

It should be noted that the derivatives blow up with respect to EPOCH and SIGN quite easily for a large signal if the model starts from a value which is far off. I don’t think this is a huge issue because I can imagine sign not actually needing to be fit. ( It should be obvious !) Similarly, it should be easy enough to start EPOCH at a nearby location. I would imagine that these should get found in a noise run first anyways.

[jeremy-baier]

Injected vs. Recovered Parameters

Parameter	Injected	Recovered	Residual	Free?
ELONG	45	45	8.168e-09	yes
ELAT	15	15	2.621e-10	yes
F0	300	300	2.842e-14	yes
F1	1e-15	1.0000008e-15	8.073e-22	yes
DM	10	9.9999957	-4.302e-06	yes
CHROMGAUSS_EPOCH_1	55000	55001.361	1.361	yes
CHROMGAUSS_LOGAMP_1	-5.30103	-5.3177291	-0.0167	yes
CHROMGAUSS_CHROMIDX_1	4	4.0643556	0.06436	yes
CHROMGAUSS_LOGSIG_1	2.4771213	2.4769894	-0.0001318	yes
CHROMGAUSS_SIGN_1	1	1	0	no

[jeremy-baier]

I think this is ready for a final review 🤠
@dlakaplan

[dlakaplan]

I finally sat down and ran a simple test of this. It seems that the model recovery happens very well with a toy model. The attached images show 1) injected chromatic Gaussian 2) intentionally slightly misspecified timing model with injected chromatic Gaussian 3).

It should be noted that the derivatives blow up with respect to EPOCH and SIGN quite easily for a large signal if the model starts from a value which is far off. I don’t think this is a huge issue because I can imagine sign not actually needing to be fit. ( It should be obvious !) Similarly, it should be easy enough to start EPOCH at a nearby location. I would imagine that these should get found in a noise run first anyways.

That's looking pretty good. But what are the 3 different plots (or presumably the last is post-fit resids, so what are the first two)?

[dlakaplan]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.
I will plan to rewrite the docs so that it matched the parameters themselves.

Great. Note that it can be frustrating to check the builds since it takes an ~hour for GitHub to do it. If you really want to dig into that I can give you instructions for doing a local build which can help iterate faster.

The docs are looking good. Nice to show an example. Can you also maybe show how that looks in the par file (print a subsection of the model object, or similar)? I am curious just because this is new, so I don't know the standard form.

[jeremy-baier]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.
I will plan to rewrite the docs so that it matched the parameters themselves.

Great. Note that it can be frustrating to check the builds since it takes an ~hour for GitHub to do it. If you really want to dig into that I can give you instructions for doing a local build which can help iterate faster.

The docs are looking good. Nice to show an example. Can you also maybe show how that looks in the par file (print a subsection of the model object, or similar)? I am curious just because this is new, so I don't know the standard form.

1. The first, top plot shows otherwise ideal residuals with a chromatic Gaussian “injected” into them. Mostly to confirm that the injection is happening. I made this by simulating residuals with a model that included the Gaussian, then removed the Gaussian from the model.
2. The 2nd plot shows an intentionally mis-specified model. ( I changed sky location, the Gaussian, and DM by a tiny amount ). This was just to mess up the timing model before refitting. This plot was probably unnecessary to share here.
3. The third plot shows the zero mean, Gaussian residuals after refitting with the Gaussian in the model.

[jeremy-baier]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.
I will plan to rewrite the docs so that it matched the parameters themselves.

Great. Note that it can be frustrating to check the builds since it takes an ~hour for GitHub to do it. If you really want to dig into that I can give you instructions for doing a local build which can help iterate faster.

The docs are looking good. Nice to show an example. Can you also maybe show how that looks in the par file (print a subsection of the model object, or similar)? I am curious just because this is new, so I don't know the standard form.

This is what the model looks like if I print the par file:

PSR                            J0000+0000
UNITS                                 TDB
DILATEFREQ                              N
DMDATA                                  N
NTOA                                    0
ELONG                  45.000000000000000 1 0.00000000000000000000
ELAT                   15.000000000000000 1 0.00000000000000000000
PMELONG                               0.0
PMELAT                                0.0
PX                                    0.0
ECL                              IERS2010
F0                                  300.0 1 0.0
F1                                  1e-15 1 0.0
PEPOCH             57000.0000000000000000
PLANET_SHAPIRO                          N
DM                                   10.0 1 0.0
CHROMGAUSS_FREF                    1400.0
CHROMGAUSS_EPOCH_1    55000.0000000000000000
CHROMGAUSS_LOGAMP_1        -5.301029995663981
CHROMGAUSS_CHROMIDX_1                       4.0
CHROMGAUSS_SIGN_1                       1.0
CHROMGAUSS_LOGSIG_1        2.4771212547196626

[dlakaplan]

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

I'll start taking a peak at this. Note that the docs render well: https://nanograv-pint--1980.org.readthedocs.build/en/1980/_autosummary/pint.models.transient_events.ChromaticGaussianEvent.html#pint.models.transient_events.ChromaticGaussianEvent but I think you could probably use more standard formatting for the paper reference (at least a URL).

Oh hey, I didn’t know you could access the docs build like that. Super handy. Thanks.
I will plan to rewrite the docs so that it matched the parameters themselves.

Great. Note that it can be frustrating to check the builds since it takes an ~hour for GitHub to do it. If you really want to dig into that I can give you instructions for doing a local build which can help iterate faster.

The docs are looking good. Nice to show an example. Can you also maybe show how that looks in the par file (print a subsection of the model object, or similar)? I am curious just because this is new, so I don't know the standard form.

This is what the model looks like if I print the par file:

PSR                            J0000+0000
UNITS                                 TDB
DILATEFREQ                              N
DMDATA                                  N
NTOA                                    0
ELONG                  45.000000000000000 1 0.00000000000000000000
ELAT                   15.000000000000000 1 0.00000000000000000000
PMELONG                               0.0
PMELAT                                0.0
PX                                    0.0
ECL                              IERS2010
F0                                  300.0 1 0.0
F1                                  1e-15 1 0.0
PEPOCH             57000.0000000000000000
PLANET_SHAPIRO                          N
DM                                   10.0 1 0.0
CHROMGAUSS_FREF                    1400.0
CHROMGAUSS_EPOCH_1    55000.0000000000000000
CHROMGAUSS_LOGAMP_1        -5.301029995663981
CHROMGAUSS_CHROMIDX_1                       4.0
CHROMGAUSS_SIGN_1                       1.0
CHROMGAUSS_LOGSIG_1        2.4771212547196626

Nice. Can you include that in the docs?

[dlakaplan]

So I think this looks largely good. Like I said I haven't verified all of the partial derivatives in detail, but the fact that fits largely converge is good. I personally am probably happy to merge it but I don't know if anybody else (@abhisrkckl, @mtlam ) wants to look

[mtlam]

Minor comments:

• In the docs, "Note a typo in Reardon et al., which lacks a negative sign in exponential" -> in "the" exponential
• Reardon et al. only has a Gaussian shape but not actually a chromatic one. I think it's a relatively straightforward construction to include this dependence as the nu^-gamma proportionality is very common, but it's worth saying in the docs in my opinion because if you go looking at that source, you won't find this chromatic Gaussian event model.
• Where does the Coles et al. paper come in here?

[jeremy-baier]

Minor comments:

• In the docs, "Note a typo in Reardon et al., which lacks a negative sign in exponential" -> in "the" exponential
• Reardon et al. only has a Gaussian shape but not actually a chromatic one. I think it's a relatively straightforward construction to include this dependence as the nu^-gamma proportionality is very common, but it's worth saying in the docs in my opinion because if you go looking at that source, you won't find this chromatic Gaussian event model.
• Where does the Coles et al. paper come in here?

The model that Reardon et al fit is dispersive:

Consequently, in addition to the Gaussian process DM
variations, we include a Gaussian-shaped DM event for
PSR J1603−7202 (Equation (7) of Goncharov et al. 2021a),
to describe its extreme scattering event (Coles et al. 2015;
Reardon & Coles 2023) and annual DM variations (Equation
(8) of Goncharov et al. 2021a) for PSR J0613−0200 (Keith
et al. 2013).

Looks like they just don’t really write the model equation up properly.

I cited Coles in the docstring because they report on extreme scattering events in pulsar timing data. They don’t use a Gaussian model there as far as I see. I will plan to update the docstring to be more clear there. Something like “See Coles et al 2015 for more details on extreme scattering events."

[jeremy-baier]

from Gonchorov et al 2021a:

Note that K is not the dm constant but just 1400 MHz so this matches the implementation here exactly I think.

[jeremy-baier]

test_chom_gauss.ipynb

I updated the docstring according to @mtlam ’s suggestion.
I have attached my test notebook in case anyone wants to try it out.

[jeremy-baier]

Checking in on this. @dlakaplan @mtlam
let me know if this needs anything else

[mtlam]

Things look fine to me.

[jeremy-baier]

Thanks to you both for taking a look at this !