BXA icon indicating copy to clipboard operation
BXA copied to clipboard

Example on how to plot different model components

Open ivvv opened this issue 3 years ago • 5 comments

I've tried to follow the docs and the API in order to visualise the different components in my model and did not succeed. I have 3 model parameters. Here is what I tried:

m = xspec.Model("constant*(bknpower + gauss)") Then I set the priors on some parameters and solved and get a nice corner plot etc. And trying to visualise the two models I'm interested in with this:

data = solver.posterior_predictions_convolved(component_names=['ignore','bknpower','gauss'],nsamples=100)

And I can only see on black curve passing through the data points, no separtion of the components. Adding at the end plt.gca().legend() only showed one legend entry for the data and nothing for the models.

I even tried to add plot_args=[{'color': 'red'},{'color:'green'},{'color':'blue'}] in the call to solver.posterior_predictions_convolved() with no success.

I suppose I'm doing something wrong.

In the same way of insufficient info in the docs, how can I get the actual total model (or per component) with the best-fit parameters? Currently it is plotted with plt.figure() but if I want to use the predicted model in subplots then it's not very convenient. I can see the dataset data['models'] with shape (100,1,1306) but I cannot make any use of it with the available info.

Can you please add some examples in the documentations?

Thanks in advance,

ivvv avatar Jun 30 '21 15:06 ivvv

Can you check what

Plot.addCommand('wdata foo') xspec.Plot('counts')

writes to the file? Those columns are what is being plotted. (expecting 6 or 4 columns, depending on Plot.background).

Is there a better way to get plot info from pyxspec?

JohannesBuchner avatar Jun 30 '21 15:06 JohannesBuchner

Yes, the two commands did write to a file foo.qdp and as far as I know, there is no better way to get this from pyXspec (I've learned this trick from Andy Beardmore).

ivvv avatar Jun 30 '21 16:06 ivvv

But I have some troubles to understand what's in the QDP file though, when used with xspec.Plot('counts').

It seems to me that the last column in the QDP file is the total model but strangely the columns for the individual components are all zeroes. Will try to understand it.

I used it for xspec.Plot('ldata')(for another problem) and the content was quite different: energy,dE,rate,rate_err,total,m1,m2,...

ivvv avatar Jun 30 '21 16:06 ivvv

Probably the code should be simplified (replacing the callback function with code duplication).

Essentially what you need is:

for k, row in enumerate(analyser.posterior):
     analyser.set_parameters(values=row, transformations=analyser.transformations)

     # do something with xspec.Plot.* and matplotlib.*

JohannesBuchner avatar Jun 30 '21 16:06 JohannesBuchner

Everything else in the code is just working around xspec weirdness

JohannesBuchner avatar Jun 30 '21 16:06 JohannesBuchner