Post by Sarah on Jul 31, 2017 15:12:48 GMT 8
Here's a couple of things I learned while trying to determine PSFs. I might add to it as I remember more...
I work with KiDS r-band data (so far), but a lot of this will probably be true for other data, too. I try to determine PSFs by finding suitable stars near my galaxy of interest and fitting them with a Moffat function using optim (downhill gradient algorithm).
1. If your model is (near) perfect, make sure to choose the normal distribution as your likelihood function (option like.func='norm' in profitSetupData, which is the default as well). Using the t-distribution might cause optim/ProFit to make fit worse than it needs to be as it desperately tries to make the residuals look like a t-distribution. That being said, you're probably better off with the t-distribution for non-perfect models (which is usually the case for more complex things, such as galaxies)
2. It is worth having a play with different algorithms, especially if interested in speed. For example, I found that the "BFGS" algorithm is a factor ~3 faster than the "L-BFGS-B" algorithm (both are options available in optim) with the same results (the difference is that "L-BFGS-B" can take limits, while "BFGS" can't; but those might not be needed if your initial guesses are good enough). Both of these are much faster than a full MCMC, which is usually not worth it for such simple cases.
3. I found that the normalisation of the cutout makes a difference in the speed and goodness of the fit. My data has a magnitude zeropoint of 0 and flux values typically around 1e-8. If I divide through by the sum of the cutout around the star (meaning I normalise the star to have a magnitude of 0 approximately; provided my background is zero and there's no other objects in the cutout), then the fit takes much longer but is also much better than if just fit the cutout without normalisation. I presume it's got something to do with numerical accuracy etc. (i.e. the absolute size of the numbers matters; and them being ~ 1 is good), although that's just a guess. I didn't find a difference for different convergence tolerance values; but I didn't look at this in much detail.
4. I tried fitting stars with just a Moffat function; or with a Moffat function and a constant background. There's a bit of a degeneracy here: If you have a background and you don't fit for it, then the Moffat function wings might be used to fit the background instead; whereas if you don't have a background and fit for it, then the wings of your star might be fit with the background instead. Obviously, both of this is unwanted. I think the best way to get around it is to adjust the size of the cutout accordingly: If not fitting a background, use only the segment of the star to fit so there's not much background at all to fit; while when fitting a background, make sure the cutout is definitely large enough to extend beyond the wings of the star (I found a minimum of ~ 20 FWHM of my star to be appropriate for my data).
5. To make my PSF estimate more robust, I use several stars around the galaxy of interest. The question then is how to combine those: The main options are fitting them individually and averaging somehow; or fitting them all together as in the Advanced PSF fit example in the Full Monty vignette. The latter may be more elegant, but can also take much longer, e.g. for one of my typical 2000x2000 pixel images with ~ 8 stars identified as being suitable, it takes about half an hour to fit them simultaneously; while fitting them individually only takes a few minutes in total. Individual fits also give an easy way to reject outliers if you're not sure your star selection is good enough - but it leaves the problem of how to meaningfully combine the remaining fits (I settled on just taking the median of the fit parameters and creating a model PSF from that; but depending on your data and time constraints, it might be nicer to use a simultaneous fit - maybe even after rejecting outliers)
6. This may be an obvious one, but always run profitSetupData again when changing any of the model parameters etc. since it does a lot of internal conversions etc. So as an example, if you run Data=profitSetupData(modellist, tofit, tolog, etc.); then decide to change the parameters you want to fit and just change Data$tofit without calling profitSetupData again, it won't work...
I work with KiDS r-band data (so far), but a lot of this will probably be true for other data, too. I try to determine PSFs by finding suitable stars near my galaxy of interest and fitting them with a Moffat function using optim (downhill gradient algorithm).
1. If your model is (near) perfect, make sure to choose the normal distribution as your likelihood function (option like.func='norm' in profitSetupData, which is the default as well). Using the t-distribution might cause optim/ProFit to make fit worse than it needs to be as it desperately tries to make the residuals look like a t-distribution. That being said, you're probably better off with the t-distribution for non-perfect models (which is usually the case for more complex things, such as galaxies)
2. It is worth having a play with different algorithms, especially if interested in speed. For example, I found that the "BFGS" algorithm is a factor ~3 faster than the "L-BFGS-B" algorithm (both are options available in optim) with the same results (the difference is that "L-BFGS-B" can take limits, while "BFGS" can't; but those might not be needed if your initial guesses are good enough). Both of these are much faster than a full MCMC, which is usually not worth it for such simple cases.
3. I found that the normalisation of the cutout makes a difference in the speed and goodness of the fit. My data has a magnitude zeropoint of 0 and flux values typically around 1e-8. If I divide through by the sum of the cutout around the star (meaning I normalise the star to have a magnitude of 0 approximately; provided my background is zero and there's no other objects in the cutout), then the fit takes much longer but is also much better than if just fit the cutout without normalisation. I presume it's got something to do with numerical accuracy etc. (i.e. the absolute size of the numbers matters; and them being ~ 1 is good), although that's just a guess. I didn't find a difference for different convergence tolerance values; but I didn't look at this in much detail.
4. I tried fitting stars with just a Moffat function; or with a Moffat function and a constant background. There's a bit of a degeneracy here: If you have a background and you don't fit for it, then the Moffat function wings might be used to fit the background instead; whereas if you don't have a background and fit for it, then the wings of your star might be fit with the background instead. Obviously, both of this is unwanted. I think the best way to get around it is to adjust the size of the cutout accordingly: If not fitting a background, use only the segment of the star to fit so there's not much background at all to fit; while when fitting a background, make sure the cutout is definitely large enough to extend beyond the wings of the star (I found a minimum of ~ 20 FWHM of my star to be appropriate for my data).
5. To make my PSF estimate more robust, I use several stars around the galaxy of interest. The question then is how to combine those: The main options are fitting them individually and averaging somehow; or fitting them all together as in the Advanced PSF fit example in the Full Monty vignette. The latter may be more elegant, but can also take much longer, e.g. for one of my typical 2000x2000 pixel images with ~ 8 stars identified as being suitable, it takes about half an hour to fit them simultaneously; while fitting them individually only takes a few minutes in total. Individual fits also give an easy way to reject outliers if you're not sure your star selection is good enough - but it leaves the problem of how to meaningfully combine the remaining fits (I settled on just taking the median of the fit parameters and creating a model PSF from that; but depending on your data and time constraints, it might be nicer to use a simultaneous fit - maybe even after rejecting outliers)
6. This may be an obvious one, but always run profitSetupData again when changing any of the model parameters etc. since it does a lot of internal conversions etc. So as an example, if you run Data=profitSetupData(modellist, tofit, tolog, etc.); then decide to change the parameters you want to fit and just change Data$tofit without calling profitSetupData again, it won't work...