Simple ProSpect Fitting using Highlander
Jessica E. Thorne
Making Mock Data
Load the libraries we need:
library(ProSpect)
library(LaplacesDemon)
library(foreach)
library(celestial)## Loading required package: RANN## Loading required package: NISTunits## Loading required package: pracma##
## Attaching package: 'pracma'## The following objects are masked from 'package:LaplacesDemon':
##
## logit, loglog, Modelibrary(magicaxis)
library(Highlander)Load some data we want to use:
data("BC03lr") #BC03 spectral library
data("Dale_NormTot") #Normalised Dale templates
data("pivwave") # Pivot/effective wavelengthsSet up using the GAMA filters, additional filters can be used from the EAZY_filters data object which comes with ProSpect
seed = 1
redshift = 0.1
agemax = 13.3e9-cosdistTravelTime(z=redshift, H0 = 67.8, OmegaM = 0.308)*1e9
filters=c('FUV_GALEX', 'NUV_GALEX', 'u_SDSS', 'g_SDSS', 'r_SDSS', 'i_SDSS', 'Z_VISTA',
'Y_VISTA', 'J_VISTA', 'H_VISTA', 'K_VISTA', 'W1_WISE' , 'W2_WISE', 'W3_WISE',
'W4_WISE', 'P100_Herschel', 'P160_Herschel', 'S250_Herschel' , 'S350_Herschel',
'S500_Herschel')
filtout=foreach(i = filters)%do%{approxfun(getfilt(i))}
temppiv=pivwave[pivwave$filter %in% filters,]Create our model galaxy We will be using a snorm_trunc SFH function and a massmap_lin metallicity with the following parameters:
inpar=c(mSFR = 0, #log-space
mpeak = 0.7, #log-space
mperiod = 0.3, #log-space
mskew = 0.3,
tau_birth = 0, #log-space
tau_screen = -0.5, #log-space
alpha_SF_birth = 1,
alpha_SF_screen = 3,
Zfinal = -2 #log-space
)Simple function to view our target SFH:
plotSFH=function(par,agemax=13.3, add=FALSE,col='black',ylim=NULL,...){
magcurve(massfunc_snorm_trunc(age=x,mSFR=10^par[1],mpeak=10^par[2],mperiod=10^par[3],
mskew=par[4], magemax=agemax),0,13.8e9,add=add,col=col,
ylim=ylim,xlab='Age (Yr)', ylab='SFR (Msol / Yr)',...)
}plotSFH(inpar)
Now we can generate our galaxies SED:
genSED=ProSpectSED(massfunc=massfunc_snorm_trunc,
mSFR=10^inpar[1],
mpeak=10^inpar[2],
mperiod=10^inpar[3],
mskew=inpar[4],
tau_birth=10^inpar[5],
tau_screen=10^inpar[6],
alpha_SF_birth=inpar[7],
alpha_SF_screen=inpar[8],
Zfinal = 10^inpar[9],
z=0.1,
Z=Zfunc_massmap_lin,
filtout=filtout,
Dale=Dale_NormTot,
speclib=BC03lr,
agemax=agemax
)At this point we create our mock photometry catalogue entry with a fractional error in flux of 0.1 assumed (roughly 0.092 mag error). The photom output of ProSpectSED is always Jansky, which is what we want for fitting (you shouldn’t fit in magnitude space!).
flux_input=data.frame(filter=temppiv$filter, pivwave=temppiv$pivwave, flux=genSED$Photom, fluxerr=genSED$Photom*0.1)
print(flux_input)## filter pivwave flux fluxerr
## 1 FUV_GALEX 1535.080 4.011918e-06 4.011918e-07
## 2 NUV_GALEX 2300.785 5.984452e-06 5.984452e-07
## 3 u_SDSS 3567.001 1.275130e-05 1.275130e-06
## 4 g_SDSS 4734.851 3.799885e-05 3.799885e-06
## 5 r_SDSS 6194.748 6.333013e-05 6.333013e-06
## 6 i_SDSS 7509.465 7.908547e-05 7.908547e-06
## 7 Z_VISTA 8833.013 9.498800e-05 9.498800e-06
## 8 Y_VISTA 10223.626 1.086040e-04 1.086040e-05
## 9 J_VISTA 12545.988 1.157941e-04 1.157941e-05
## 10 H_VISTA 16476.995 1.312316e-04 1.312316e-05
## 11 K_VISTA 21549.195 1.208860e-04 1.208860e-05
## 12 W1_WISE 33897.093 7.034516e-05 7.034516e-06
## 13 W2_WISE 46406.422 4.483248e-05 4.483248e-06
## 14 W3_WISE 125700.270 2.740481e-04 2.740481e-05
## 15 W4_WISE 223141.602 6.297586e-04 6.297586e-05
## 16 P100_Herschel 983273.314 7.548473e-03 7.548473e-04
## 17 P160_Herschel 1546806.337 9.272729e-03 9.272729e-04
## 18 S250_Herschel 2482806.324 6.245763e-03 6.245763e-04
## 19 S350_Herschel 3484355.144 3.030799e-03 3.030799e-04
## 20 S500_Herschel 5004155.150 1.060366e-03 1.060366e-04Fitting a Galaxy
In principle if you have a real galaxy to fit and it is already in the required format of flux_input, you can start from this point.
To speed things up we will also pre-compute the Luminosity distance to our source (since this will be the same for every iteration, and will save some computation time):
LumDist_Mpc = cosdistLumDist(z=0.1, H0 = 67.8, OmegaM = 0.308)Data=list(flux=flux_input,
arglist=list(z=0.1,
massfunc=massfunc_snorm_trunc,
agemax=agemax,
Zagemax = (agemax/1e9),
magemax=(agemax/1e9),
Z=Zfunc_massmap_lin,
LumDist_Mpc=LumDist_Mpc),
speclib=BC03lr,
Dale=Dale_NormTot,
filtout=filtout,
Dale_M2L_func=Dale_M2L_func, # required to get monitored dust masses
SFH=SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
parm.names=c('mSFR','mpeak','mperiod','mskew','tau_birth','tau_screen', 'alpha_SF_birth','alpha_SF_screen', 'Zfinal'), # which parameters to fit for
mon.names=c("LP","masstot","dustmass.birth", "dustmass.screen", "dustmass.total", "dustlum.birth", "dustlum.screen", "dustlum.total", "SFRburst", paste("flux.",filters,sep='')), # some of the things you may wish to monitor in each step of the chain
logged=c(T,T,T,F,T,T,F,F,T), # fit parameters in logged or linear space
intervals=list(lo=c(-4,-2,-1,-0.5,-2.5,-2.5,0,0,-4),
hi=c(3,1,1,1,1.5,1,4,4,-1.3)), # fitting range for parameters
fit = 'LD', # specifies the way in which the SED should be fitted ('LD', 'optim', 'CMA', or 'check')
mon.names=c('LP','masstot','SFRburst',paste('flux.',flux_input$filter,sep='')),
N=length(filters), # number of observed filters
like='st', # Using a student-t likelihood
verbose=FALSE
)Important things:
- agemax is in years whereas magemax and Zagemax are in Gyrs - make sure you don’t use 13.8 for agemax otherwise your SED will be forced to fit all the SFH into the last 13.8 years.
- for the monitored parameters to actually be monitored you need to be using fit = ‘LD’ in the Data structure and likefunctype = ‘LD’ in the call to Highlander
Using Highlander for the fits:
startpoint = (Data$intervals$lo+Data$intervals$hi)/2
testHigh = Highlander(startpoint, Data,
ProSpectSEDlike, Niters=c(200,200), NfinalMCMC = 200,
lower=Data$intervals$lo, upper=Data$intervals$hi,
seed=seed, optim_iters = 2, likefunctype = 'LD')## Iteration 1## CMA 1: -27.687 -0.124 0.768 0.712 -0.279 1.384 -0.747 2.706 1.956 -1.906##
## Laplace's Demon was called on Fri Aug 13 12:13:35 2021
##
## Performing initial checks...
## Algorithm: Componentwise Hit-And-Run Metropolis
##
## Laplace's Demon is beginning to update...
## Iteration: 100, Proposal: Componentwise, LP: -31.8
## Iteration: 200, Proposal: Componentwise, LP: -27.9
##
## Assessing Stationarity
## Assessing Thinning and ESS
## Creating Summaries
## Creating Output
##
## Laplace's Demon has finished.## LD Mode 1: -26.461 -0.136 0.752 0.652 0.016 0.593 -0.62 2.459 0 -1.779## Iteration 2## CMA 2: -23.502 -0.121 0.479 0.492 -0.389 0.379 -0.545 2.603 0.114 -1.992##
## Laplace's Demon was called on Fri Aug 13 12:14:23 2021
##
## Performing initial checks...
## Algorithm: Componentwise Hit-And-Run Metropolis
##
## Laplace's Demon is beginning to update...
## Iteration: 100, Proposal: Componentwise, LP: -27.5
## Iteration: 200, Proposal: Componentwise, LP: -28.6
##
## Assessing Stationarity
## Assessing Thinning and ESS
## Creating Summaries
## Creating Output
##
## Laplace's Demon has finished.## LD Mode 2: -23.5 -0.121 0.479 0.492 -0.389 0.355 -0.545 2.603 0.164 -1.992Data$fit = 'check'
bestfit=ProSpectSEDlike(testHigh$par, Data=Data)Checking the output
plot(bestfit)