Bug?: BasicPSFPhotometry PSF flux uncertainty shrinks with increasing fitsize

[emirkmo]

This was previously reported, but the fitsize parameter of BasicPSFPhotometry has a shrinking uncertainty with increasing fitsize. I'm attaching a figure which shows this.

To be clear what we are doing is getting a star cutout (region of image with a few x FWHM, centered around the star [we often have sparse fields]), using our EPSF (built previously), and running BasicPSFPhotometry on the input to get a PSF flux and flux uncertainty.

The PSF flux uncertainty (error bars) shrinks with increasing fitsize even though the determined flux reaches an absolute peak value pretty quickly at around ~2x the image FWHM. This is also easily visible in the full calibrated instrumental magnitude plot attached at the end. The most likely reason for this decreasing uncertainty is some sort of overfitting, or overconfidence of the fit to the relatively smooth background. This creates an issue where the error bars are smaller than they should be. Basically, the background, instead of adding uncertainty via Poisson noise, takes away from the error budget as more of it is included in the fit.

Is there a smart way of getting around this while continuing to use the flux uncertainty of BasicPSFPhotometry, without finding the "elbow" manually for each PSF fit where the uncertainty begins to shrink while the flux does not change?

The uncertainty comes from the fit_info from the fitter. The basic call signature is:

photometry_obj = BasicPSFPhotometry(group_maker=DAOGroup(fwhm), bkg_estimator=MedianBackground(), psf_model=epsf,
                                        fitter=fitting.LevMarLSQFitter(), fitshape=phot.star_size, aperture_radius=fwhm)
    psfphot_tbl = photometry_obj(image=image.subclean, init_guesses=Table(clean_references.xy, names=['x_0', 'y_0']))

Basically the fitter is fitting.LevMarLSQFitter(), fitshape was by default the full cutout size, but has now been made a multiple of the fwhm (previous determined from the image and is in units of image pixels), and DAOGroup helps with stellar crowding, and bkg_estimator is MedianBackground BUT the image has already been background subtracted previously using another estimator such as BackGround2D.SExtractorBackground with sigma clipping.

Finally, initial guesses for the stellar positions (PSF centers) are passed in, although it is not forced, meaning the PSF can move around.

I wonder if the bkg_estimator is doing something wonky due to its smoothness? Or, I suspect that the uncertainty estimate is actually poorly defined because we are interested in not how well the PSF fits the relatively smooth background at multiple FWHM but in how it fits the majority of the flux at 2-3 FWHM from the center, and that is the uncertainty we care about.

Previous issue see: #1130

[larrybradley]

I'd need to investigate this further. The PSF flux uncertainties come from the astropy fitter class which in turn come from the scipy optimizers, e.g., scipy.optimize.leastsq for the LevMarLSQFitter fitter.

[emirkmo]

@larrybradley I think that Astropy 5 may have solved one aspect of this issue, where nans/masked values were involved. But the behavior still remains. Just an update.

I think it may be a conceptual issue of how we define the error and approximating it via gaussian statistics also in the background region. Perhaps PSF fitter could take in an additional parameter to find this elbow and re-calculate the error using it, thereby deriving a scaling factor for the uncertainty based on the fitsize but really only taking into account regions with significant background flux. We are actually already doing this with our flows photometry pipeline, see below.

Alternatively, if we could re-compute the correlation matrix while masking out all regions identified as background, then I suppose the error should be more reflective of the flux falling on relevant pixels. Would need to add a threshold parameter for this, such as 10% over median or whatever background identifier used. It's probably a more robust way of doing it for PSFs that deviate from a simple single Gaussian.

    def rescale_uncertainty(self, psfphot_tbl: Table, dynamic: bool = True, 
                            static_fwhm: float = 2.5, epsilon_mag: float = 0.004,
                            ensure_greater: bool = True):
        """
        Rescale the uncertainty of the PSF photometry using a variable fitsize.

        Parameters
        ----------
        psfphot_tbl : Table
            Table of PSF photometry.
        dynamic : bool
            Dynamically decide FWHM multiple for rescaling.
        static_fwhm : float
            FWHM multiple to use incase dynamic fails or don't want to use it. Default 2.5 determined empirically.
        epsilon_mag : float
            Small magnitude change within which new and old uncertainties are considered the same. 
            Should be smaller than ~1/2 the expected uncertainty.
        """
        # Rescale psf errors from fit iteratively
        fit_shapes = self.photometry.get_fit_shapes(self.fwhm, self.psf_builder.star_size)
        fit_shape = int(static_fwhm * self.fwhm)
        fit_shape = fit_shape if fit_shape % 2 == 1 else fit_shape + 1
        if dynamic and len(fit_shapes) > 1:
            _phot_tables_dict = {}
            for fitshape in fit_shapes:
                self.photometry.create_photometry_object(
                    fwhm=self.fwhm, psf_model=self.psf_builder.epsf, fitsize=fitshape)
                if self.diff_im_exists:
                    _table = self.diff_psf_phot()
                _table = self.raw_psf_phot(self.init_guesses.init_guess_target)
                if "flux_unc" in _table.colnames:
                    _phot_tables_dict[fitshape] = _table

            if len(_phot_tables_dict) == 0:
                logger.warning("No PSF errors found for dynamic rescaling, trying static.")
                dynamic = False
            else:
                # Find the fit shape elbow:
                flux = psfphot_tbl[0]['flux_fit']
                flux_err = psfphot_tbl[0]['flux_unc']
                exptime = self.image.exptime
                dynamic_fit_shape, new_err = self.photometry.rescale_flux_error(_phot_tables_dict, flux, flux_err,
                                                                                exptime)
                fit_shape = dynamic_fit_shape if dynamic_fit_shape != 0 else fit_shape

        # Recalculate all reference uncertainties using new fitsize:
        logger.info(f"Recalculating all reference uncertainties using new fitsize:"
                    f" {fit_shape} pixels, ({fit_shape/self.fwhm if dynamic else static_fwhm :.2} * FWHM).")
        psfphot_tbl_rescaled = self.psfphot(fit_shape)
        if psfphot_tbl['flux_unc'][0] > psfphot_tbl_rescaled['flux_unc'][0] + epsilon_mag and ensure_greater:
            logger.info("Recalculated uncertainties were smaller than original and ``ensure_greater`` was True:"
                        "Not using rescaled uncertainties for the SN.")
            psfphot_tbl['flux_unc'][1:] = psfphot_tbl_rescaled['flux_unc'][1:]
            return psfphot_tbl

        psfphot_tbl['flux_unc'] = psfphot_tbl_rescaled['flux_unc']
        return psfphot_tbl

https://github.com/SNflows/flows/blob/c7e98d25e6b5c84248c8cb1ae9dad7013837178a/flows/photometry.py#L319-L375

[larrybradley]

@emirkmo Error arrays can now be input into PSFPhotometry. The error arrays are used to weight the fits, and when used with astropy >= 5.3, they are propagated to the parameter uncertainties. Can you please try with the new PSFPhotometry class and see if it resolves this issue?

[larrybradley]

Please reopen if you still have issues. Thanks.