Setup the fitting data structures
For our initial model we treat component 1 as the putative bulge and component 2 as the putative disk. We are going to attempt a fit where the disk is forced to have nser=1 and the bulge has an axial ratio of 1.
modellist=list(
sersic=list(
xcen= c(dim(image)[1]/2, dim(image)[1]/2),
ycen= c(dim(image)[2]/2, dim(image)[2]/2),
mag= c(ExampleInit$sersic.mag1[useloc], ExampleInit$sersic.mag2[useloc]),
re= c(ExampleInit$sersic.re1[useloc], ExampleInit$sersic.re2[useloc])*
if(datasource=='KiDS'){1}else{0.2/0.339},
nser= c(ExampleInit$sersic.nser1[useloc], 1), #Disk is initially nser=1
ang= c(ExampleInit$sersic.ang2[useloc], ExampleInit$sersic.ang2[useloc]),
axrat= c(1, ExampleInit$sersic.axrat2[useloc]), #Bulge is initially axrat=1
box=c(0, 0)
)
)
modellist## $sersic
## $sersic$xcen
## [1] 84.5 84.5
##
## $sersic$ycen
## [1] 93.5 93.5
##
## $sersic$mag
## [1] 16.83217 16.83217
##
## $sersic$re
## [1] 7.0574 14.1148
##
## $sersic$nser
## [1] 4.3776 1.0000
##
## $sersic$ang
## [1] 140.8191 140.8191
##
## $sersic$axrat
## [1] 1.0000 0.4891
##
## $sersic$box
## [1] 0 0The pure model (no PSF):
The original image:
The convolved model (with PSF):
Next we define our list of what we want to fit (where TRUE means we will fit it later):
tofit=list(
sersic=list(
xcen= c(TRUE,NA), #We fit for xcen and tie the two togther
ycen= c(TRUE,NA), #We fit for ycen and tie the two togther
mag= c(TRUE,TRUE), #Fit for both
re= c(TRUE,TRUE), #Fit for both
nser= c(TRUE,FALSE), #Fit for bulge
ang= c(FALSE,TRUE), #Fit for disk
axrat= c(FALSE,TRUE), #Fit for disk
box= c(FALSE,FALSE) #Fit for neither
)
)Now we define what parameters should be fitted in log space:
tolog=list(
sersic=list(
xcen= c(FALSE,FALSE),
ycen= c(FALSE,FALSE),
mag= c(FALSE,FALSE),
re= c(TRUE,TRUE), #re is best fit in log space
nser= c(TRUE,TRUE), #nser is best fit in log space
ang= c(FALSE,FALSE),
axrat= c(TRUE,TRUE), #axrat is best fit in log space
box= c(FALSE,FALSE)
)
)Now we specify the priors object. The priors object is really a function that takes modellist as an input and returns an addative log-likelihood sum that increments the total log-likelihood for the current realisation of the model.
The elements you access within the modellist structure must really exist or ProFit will stop with an error. The modellist used internally will be the linear (un-logged version), so logging must be made within the function if, e.g., log-normal distributions are required. You do not have to apply priors to each parameter- notice below a few are left out.
priors=function(new,init,sigmas=c(2,2,2,2,5,5,1,1,1,1,30,30,0.3,0.3)){
LL=sum(
dnorm(new$sersic$xcen,init$sersic$xcen,sigmas[1:2],log=TRUE),
dnorm(new$sersic$ycen,init$sersic$ycen,sigmas[3:4],log=TRUE),
dnorm(new$sersic$mag,init$sersic$mag,sigmas[5:6],log=TRUE),
dnorm(log10(new$sersic$re),log10(init$sersic$re),sigmas[7:8],log=TRUE),
dnorm(log10(new$sersic$nser),log10(init$sersic$nser),sigmas[9:10],log=TRUE),
dnorm(log10(new$sersic$axrat),log10(init$sersic$axrat),sigmas[13:14],log=TRUE)
)
return(LL)
}
print(priors(modellist,modellist))## [1] -14.61078Notice that the values for sigmas were defined along with the function definition. In practice, you may want to define them separately as such:
sigmavec=c(2,2,5,1,1,30,0.3,0.2)
priors=function(new,init,sigmas=sigmavec){
return(sum(2*dnorm(log10(new$sersic$re),log10(init$sersic$re),sigmas[4],log=TRUE)))
}
print(priors(modellist,modellist))## [1] -3.675754Unfortunately, R does not store the provided values of sigmas, but simply the name of a variable sigma, so if you happen to redefine sigma later on, the prior function will change (or fail entirely if you delete sigma or save and then load a fit):
## [1] -0.9031654## function(new,init,sigmas=sigmavec){
## return(sum(2*dnorm(log10(new$sersic$re),log10(init$sersic$re),sigmas[4],log=TRUE)))
## }
## <bytecode: 0x7fe374c2c370>Because of this issue, we provide a utility function to create a valid prior function, saving the values of the prior distributions:
sigmas=c(2,2,5,1,1,30,0.3,Inf)
sigmas=list(
sersic=list(
xcen= numeric(2)+sigmas[1],
ycen= numeric(2)+sigmas[2],
mag= numeric(2)+sigmas[3],
re= numeric(2)+sigmas[4],
nser= numeric(2)+sigmas[5],
ang= numeric(2)+sigmas[6],
axrat= numeric(2)+sigmas[7],
box= numeric(2)+sigmas[8]
)
)
priors=profitMakePriors(modellist, sigmas, tolog, allowflat=TRUE)
print(priors(modellist,modellist))## [1] -23.25105## function (new, init, sigmas = c(sersic.xcen1 = 2, sersic.xcen2 = 2,
## sersic.ycen1 = 2, sersic.ycen2 = 2, sersic.mag1 = 5, sersic.mag2 = 5,
## sersic.re1 = 1, sersic.re2 = 1, sersic.nser1 = 1, sersic.nser2 = 1,
## sersic.ang1 = 30, sersic.ang2 = 30, sersic.axrat1 = 0.3, sersic.axrat2 = 0.3,
## sersic.box1 = Inf, sersic.box2 = Inf), tolog = c(sersic.xcen1 = FALSE,
## sersic.xcen2 = FALSE, sersic.ycen1 = FALSE, sersic.ycen2 = FALSE,
## sersic.mag1 = FALSE, sersic.mag2 = FALSE, sersic.re1 = TRUE,
## sersic.re2 = TRUE, sersic.nser1 = TRUE, sersic.nser2 = TRUE,
## sersic.ang1 = FALSE, sersic.ang2 = FALSE, sersic.axrat1 = TRUE,
## sersic.axrat2 = TRUE, sersic.box1 = FALSE, sersic.box2 = FALSE
## ), tofit = NULL, means = unlist(init), allowflat = TRUE)
## {
## LL = 0
## parms = unlist(new)
## if (!is.null(tofit))
## ps = which(tofit)
## else ps = 1:length(parms)
## for (p in ps) {
## if (!(allowflat && (sigmas[p] == Inf))) {
## parm = parms[[p]]
## mean = means[[p]]
## if (tolog[p]) {
## parm = log10(parm)
## mean = log10(mean)
## }
## LL = LL + dnorm(parm, mean, sigmas[p], log = TRUE)
## }
## }
## return = LL
## }
## <bytecode: 0x7fe3654d0a30>The hard intervals should also be specified in linear space:
intervals=list(
sersic=list(
xcen=list(lim=c(0,300),lim=c(0,300)),
ycen=list(lim=c(0,300),lim=c(0,300)),
mag=list(lim=c(10,30),lim=c(10,30)),
re=list(lim=c(1,100),lim=c(1,100)),
nser=list(lim=c(0.5,20),lim=c(0.5,20)),
ang=list(lim=c(-180,360),lim=c(-180,360)),
axrat=list(lim=c(0.1,1),lim=c(0.1,1)),
box=list(lim=c(-1,1),lim=c(-1,1))
)
)It is not a requirement to have a constraints object that do more complex manipulations of the modellist, but it is useful to stop undesirable (perhaps unphysical) exploration.
The constraint object is really a function that takes modellist as an input and returns a modifed (but identical in skeleton structure) modellist output. For this reason it is generally best to make manipulations on the modellist directly rather than attempting to unlist (collapse) it and rebuild it etc (this can be done, but it is prone to coding error).
The elements you access within the modellist structure must really exist or ProFit will stop with an error. All constraints manipulations happen internally to the model in linear (un-logged) form.
Setup the data structure we need for optimisation, taking a few seconds to find the optimal convolution method:
Data=profitSetupData(image=image, sigma=sigma, segim=segim, psf=psf,
modellist=modellist, tofit=tofit, tolog=tolog, priors=priors,
intervals=intervals,magzero=0, algo.func='optim', like.func="t",
verbose=TRUE)This produces a fairly complex R object, but with all the bits we need for fitting, e.g. (notice the tolog parameteres are now logged):
## sersic.xcen1 sersic.ycen1 sersic.mag1 sersic.mag2 sersic.re1
## 84.5000000 93.5000000 16.8321750 16.8321750 0.8486447
## sersic.re2 sersic.nser1 sersic.ang2 sersic.axrat2
## 1.1496747 0.6412361 140.8191000 -0.3106023These are the parameters we wish to fit for, and we take the initial guesses from the model list we provided before.
We can test how things currently look (we get an output because we set verbose=TRUE earlier):
## [1] "Summary of used sample:"
## Min. 1st Qu. Median Mean 3rd Qu. Max.
## -76.6588 -0.4123 0.4049 0.2498 1.2617 23.7657
## [1] "sd / 1-sig / 2-sig range:"
## [1] 3.558640 1.280227 4.606598
## [1] "Using 16777 out of 16777"
## sersic.xcen1 sersic.ycen1 sersic.mag1 sersic.mag2 sersic.re1
## 8.4500e+01 9.3500e+01 1.6832e+01 1.6832e+01 8.4864e-01
## sersic.re2 sersic.nser1 sersic.ang2 sersic.axrat2
## 1.1497e+00 6.4124e-01 1.4082e+02 -3.1060e-01 -3.3484e+04