Linear Timing Model Coefficients

[Hazboun6]

As part of work on the nonlinear timing model functionality @ark0015 and I have been working on a version of the linear timing model that allows the user to choose which timing parameters will be modeled with GPs in the linear approximation. This is done by slicing the design matrix according to the list of parameter names. The hope is that this will allow for a bit more freedom in the modeling and allow the parameters of interest in the nonlinear timing model to explore the parameter space more effectively. The functionality currently works when we are not concerned with keeping the linear timing model perturbations, but breaks when we set the coefficients=True flag in TimingModel. We'd like to be able to keep the coefficients in order to investigate if other parameters are covariant, or otherwise change significantly when doing the nonlinear timing model fit.

We get the following error:
AttributeError: No sampler was provided for this Parameter.
This appears to stem from the lack of all of the parameter attributes when trying to use the sample function.

We found the following test for the GPCoefficient functionality, but trying to use the utils.get_coefficients(pta, ps) function and putting the results into the likelihood calculation return the same error.

There is a "ToDo" listed in the code by @vallis.

@vallis and @stevertaylor I'm wondering if this functionality exists in one of your branches for the cross validation study, or perhaps it wasn't needed? How hard would it be to implement the changes needed to keep the linear timing model coeefficients?

[vallis]

Sorry for missing this… the easiest way to exclude parameters from the marginalization and “coefficients” mechanism is to turn them off with libstempo, e.g., psr[‘RAJ’].fit = False before the design matrix is saved into the pulsar object/enterprise. If parameters are not modeled linearly, then it does not make sense to include them in “coefficients”. This is just a quick answer without looking at code, but perhaps enough to continue the discussion? Michele

…

On Sep 10, 2020, 2:20 PM -0700, Jeff Hazboun ***@***.***>, wrote: As part of work on the nonlinear timing model functionality @ark0015 and I have been working on a version of the linear timing model that allows the user to choose which timing parameters will be modeled with GPs in the linear approximation. This is done by slicing the design matrix according to the list of parameter names. The hope is that this will allow for a bit more freedom in the modeling and allow the parameters of interest in the nonlinear timing model to explore the parameter space more effectively. The functionality currently works when we are not concerned with keeping the linear timing model perturbations, but breaks when we set the coefficients=True flag in TimingModel. We'd like to be able to keep the coefficients in order to investigate if other parameters are covariant, or otherwise change significantly when doing the nonlinear timing model fit. We get the following error: AttributeError: No sampler was provided for this Parameter. This appears to stem from the lack of all of the parameter attributes when trying to use the sample function. We found the following test for the GPCoefficient functionality, but trying to use the utils.get_coefficients(pta, ps) function and putting the results into the likelihood calculation return the same error. There is a "ToDo" listed in the code by @vallis. @vallis and @stevertaylor I'm wondering if this functionality exists in one of your branches for the cross validation study, or perhaps it wasn't needed? How hard would it be to implement the changes needed to keep the linear timing model coeefficients? — You are receiving this because you were mentioned. Reply to this email directly, view it on GitHub, or unsubscribe.

[Hazboun6]

Hi @vallis, Thanks for the response. It's good to know that we can turn off the fit in the libstempo object. I assume we would then need to regenerate the design matrix in the enterprise Pulsar object then, is that correct?

The more important issue that we are having is that if we try and use the coefficient functionality, without even changing anything else like which parameters are fit, we still can't use pta.get_lnlikelihood. We can't even construct an initial sample.

When you have used this functionality in the past do you sample with these coefficients turned on? Are they written to the chains? As stated in the note I linked above, there is no sample method for the GPCoefficients, so when we try to calculate an initial sample we get an exception. Also if we try and use the get_coefficients method we also get an error:
AttributeError: 'GPCoefficients' object has no attribute 'value'

[Hazboun6]

@vallis This is the issue that I was discussing on the DWG telecon yesterday. The issue is basically that it would be very useful in @ark0015 's work to be able to recover the timing model parameter changes from the linearized timing model in order to compare to the nonlinear timing model parameters that he is sampling over explicitly. It seems like there are other instances when this would be useful functionality for tracking down some of these instances where unconverged timing models are used.

I am not sure this is functionality needed for 15 year analyses, but as I said it would definitely be helpful for the nonlinear timing comparisons. There are also use cases in the advanced noise modeling where it would be useful to know when the timing parameters are changing substantially when a chromatic GP is used.

[vallis]

Hello @Hazboun6 and @ark0015, sorry for not answering back in September. The modus operandi for get_coefficients is to feed it hyperparameter sets with a PTA with coefficients=False; the function then returns a dictionary of parameters that include the coefficient vectors, and that can be used with another PTA that has coefficients=True. (I used this for posterior likelihoods, so I am not sure (I don't remember) whether it's possible to MCMC-sample the PTA with coefficients=True.)

If you want to see the implicit value of timing model parameters with the marginalized statistic, you can almost do it with get_coefficients. You'd want to have gp_signals.TimingModel(use_svd=False, normed=False) in the PTA so the design-matrix vectors don't get mixed up and the coefficient units are consistent. Then get_coefficients will give you values for parameters such as B1855+09_linear_timing_model_coefficients; those values will be sampled from the conditional distribution of the coefficients given the hyperparameters.

The conditional is jointly normal, and you can get multiple samples with the option n. If you want the conditional mean, you can with a small code change, see my branch https://github.com/vallis/enterprise/tree/getcoefficients. With a little more work I could give you the covariance if you need it.

Let me know if this makes sense and if gives sensible results.

[ark0015]

Thanks for looking into this @vallis, it's been awhile since I've looked at what we tried before, but just running what you suggest gives the linear_timing_model_coefficients param, but still gives an error. This time it is giving ValueError: could not broadcast input array from shape (0) into shape (87), which I think is because the initial sample may not provide draws for all the timing model parameters.

[vallis]

What did you run exactly, @ark0015? If you used my new branch, it's possible I botched the tweak (the equivalent of missing the easy putt).

But if you get a vector of linear_timing_model_coefficients, isn't that what you want?

[ark0015]

@vallis, I setup a single pulsar pta with WN and with gp_signals.TimingModel(use_svd=False, normed=False, coefficients=True), then I was trying to mimic this test for the GPCoefficients and when doing psc = utils.get_coefficients(pta, x0_dict), I get enterprise/signals/utils.py, line 88, in get_coefficients Sigma = TNT + (np.diag(phiinv) if phiinv.ndim == 1 else phiinv) AttributeError: 'NoneType' object has no attribute 'ndim'. I'm not really sure why get_phiinv is returning None, but could the issue be in there?

Ideally we would be able to get the linear timing values back at each sample to determine whether they are moving into nonlinear territory, is that feasible with this kind of setup?

[vallis]

get_coefficients needs to be called on a PTA with coefficients=False (the default). The idea is that we sample the marginalized posterior, and then use the conditional posterior of the coefficients (given the data and hyperparameters) to sample their values. With the patch in my branch, get_coefficients can return their mean instead of a random value from their normal distribution.