Total flux error propagation in joint (Sersic) models

[entaylor]

Describe the bug

I am fitting a single Sersic profile to 6 images of the same galaxy, following the example of the 'Joint Modelling' tutorial. More specifically, i create Target_Image and AstroPhot_Model objects for each image, with the model parameters for all the other models locked to the values of the first like this:

                    model0 = model
                else :
                    for p in joint_params :
                        model[p].value = model0[p]
                model_list.append( model )

I then combine these into a group model with the list of models and the targets packaged in a Target_Image_List object.

The fitting process all seems to proceed correctly (or as expected!) After convergence and updating uncertainties, i can get the model parameters and their uncertainties, which all look reasonable:

q: 0.870119221222134 +- 0.0007164113555186324 [b/a], limits: (0.0, 1.0)
PA: 1.8315863639907133 +- 0.0029465403669712846 [radians], limits: (0.0, 3.141592653589793), cyclic
n: 2.5392492286223116 +- 0.002599553314399845 [none], limits: (0.36, 8.0)
Re: 83.43737384049845 +- 0.14919776412953098 [arcsec], limits: (0.0, None)
Ie: 0.32539842725159474 +- 0.0012260773462222482 [log10(flux/arcsec^2)]
Ie: 0.3217348811739831 +- 0.0012268533564793176 [log10(flux/arcsec^2)]
Ie: 0.32202030241809004 +- 0.0012267924450888779 [log10(flux/arcsec^2)]
Ie: 0.3163010787086617 +- 0.0012280277887813905 [log10(flux/arcsec^2)]
Ie: 0.32448974030480787 +- 0.0012262686501713599 [log10(flux/arcsec^2)]
Ie: 0.32516677247767256 +- 0.0012261260428432467 [log10(flux/arcsec^2)]

This is a very big and bright galaxy, so Re is large but everything is very well constrained. So far so good.

Expected behavior

The problem arises when i try to use the total_flux and total_flux_uncertainty methods to extract the total flux and associated uncertainty. Trying to call these methods attached to the group model fails with a type error, which i think makes sense: presumably, the total flux of the group model would add together the total flux in each image, which is not what i want. Calling the methods attached to each of the models that goes into the group works, and gives sensible results for the total fluxes (hooray!) but not for the uncertainties. If i:

            sersic_flux_unc = np.array([ model.total_flux_uncertainty().item() for model in model_full ])
            print( sersic_flux )
            print( sersic_flux_unc )

the results are:

[233757.87861623 231794.27818483 231946.66500726 228912.18394143
 233269.29154434 233633.22428349]
[1.08820865e+03 1.43611482e+08 1.53607006e+09 7.43031667e+07
 2.19916433e+08 1.14584884e+08]

For testing, i am using 6 images in the same band, so the fluxes should all be the same (and they are) and the uncertainties should all be comparable (which they are not). You can see this, too, comparing the errors in Ie to the reported errors in total flux.

Given that i have plausible errors in the Ies and all the other parameters, i can do the propagation of uncertainties to get my own values for the uncertainties in total flux, but ... hopefully it's also easy to fix this in the code?

You're amazing!

[ConnorStoneAstro]

Hi @entaylor this is very interesting. Yes, it looks like everything is working up to the point of the total flux uncertainties and even then it looks like the first one is ok (total flux is 2.3e5 and uncertainty is 1e3) and the others are crashing. Could you also plot the target, model, and residual images? Do they look substantially different from each other?

Also just to check, I assume you are using the latest version? I made an update to how the uncertainty is calculated and the formula is quite simple (but that has never stopped me before from introducing a bug, haha).

[ConnorStoneAstro]

Along the same lines, are you providing a variance image for each target? They need that to be accurate in order to compute uncertainties correctly.

[entaylor]

This was with 0.16.12; i've just upgraded to 0.16.13 and ... now that works as expected!! So thank you! Does that make sense to you that something relevant may have changed in that upgrade?

While i'm here: yep, i'm feeding good variance maps (and get the same results if i instead pass inverse-variance weights). Just for fun, my variance maps now include shot noise, which i have found to be very important in the joint PSF fits i've posted about previously. I took me a while to understand this, but ... yeah, accounting for shot noise really makes a difference to joint fitting, and it'd be fun to talk over some time if you have the interest and opportunity. : )

Then, because i have them, here are the pretty joint fits to six independent images:

[ConnorStoneAstro]

Hi @entaylor that's great! Yes v0.16.13 had an update to how uncertainties are calculated, I changed it to a very general formula rather than relying on custom, manual derivatives.

Oh I totally agree that better variance maps (e.g. including shot noise like your case) is critical for multi image fitting. The variance tells the likelihood about the relative weights of the images so its important that those are correct for the multiple images to interact properly. I'd be happy to chat more, but I can say that the ideal case is when a variance or weight map is reported alongside the data. I've been working with some HST data lately and they are really good for this.

Those are very pretty images!