Making Mock Data
Load the libraries we will 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, Mode
Load some data we want to use:
data("BC03lr") #BC03 spectral library
data("Dale_NormTot") #Normalised Dale templates
data("pivwave") # Pivot/effective wavelengths
set.seed(1)
redshift = 0.1
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,]
agemax = 13.3e9-cosdistTravelTime(z=redshift, H0 = 67.8, OmegaM = 0.308)*1e9
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
)
Simple function to view our target SFH:
plotSFH=function(par, magemax=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=magemax),0,13.8e9,add=add,col=col,
ylim=ylim,xlab='Age (Yr)', ylab='SFR (Msol / Yr)',...)
}
Let’s take a look:
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],
z=0.1,
Z=Zfunc_massmap_lin,
filtout=filtout,
Dale=Dale_NormTot,
speclib=BC03lr,
agemax=agemax,
magemax=agemax/1e9
)
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 3.295781e-06 3.295781e-07
## 2 NUV_GALEX 2300.785 5.120805e-06 5.120805e-07
## 3 u_SDSS 3567.001 1.084647e-05 1.084647e-06
## 4 g_SDSS 4734.851 3.252644e-05 3.252644e-06
## 5 r_SDSS 6194.748 5.551887e-05 5.551887e-06
## 6 i_SDSS 7509.465 7.030547e-05 7.030547e-06
## 7 Z_VISTA 8833.013 8.615903e-05 8.615903e-06
## 8 Y_VISTA 10223.626 1.013659e-04 1.013659e-05
## 9 J_VISTA 12545.988 1.115099e-04 1.115099e-05
## 10 H_VISTA 16476.995 1.290294e-04 1.290294e-05
## 11 K_VISTA 21549.195 1.226153e-04 1.226153e-05
## 12 W1_WISE 33897.093 7.085921e-05 7.085921e-06
## 13 W2_WISE 46406.422 4.586332e-05 4.586332e-06
## 14 W3_WISE 125700.270 2.447705e-04 2.447705e-05
## 15 W4_WISE 223141.602 5.847274e-04 5.847274e-05
## 16 P100_Herschel 983273.314 6.831858e-03 6.831858e-04
## 17 P160_Herschel 1546806.337 8.199005e-03 8.199005e-04
## 18 S250_Herschel 2482806.324 5.474332e-03 5.474332e-04
## 19 S350_Herschel 3484355.144 2.651837e-03 2.651837e-04
## 20 S500_Herschel 5004155.150 9.275196e-04 9.275196e-05
Fitting 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, magemax=agemax/1e9,
Z=Zfunc_massmap_lin, LumDist_Mpc=LumDist_Mpc),
speclib=BC03lr,
Dale=Dale_NormTot,
filtout=filtout,
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'), # which parameters to fit for
logged=c(T,T,T,F,T,T,F,F), # fit parameters in logged or linear space
intervals=list(lo=c(-4,-2,-1,-0.5,-2.5,-2.5,0.0625,0.0625), hi=c(3,1,1,1,1,1,4,4)), # 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='norm',
verbose=FALSE
)
Now we can run an MCMC fit using LaplacesDemon (depending on machine, this might take 10-20 minutes):
set.seed(1)
LDout = LaplacesDemon(Model=ProSpectSEDlike, Data=Data, Initial.Values=inpar,
control=list(abstol=0.1), Iterations=1e4, Algorithm='CHARM', Thinning=1)
Above we cheated by starting on the right values for the fit. In reality we would have to use an optimiser to get a reasonable starting position. R comes with many, but optim is a good simple start:
set.seed(1)
Data$fit = 'optim'
optout = optim(par=inpar, fn=ProSpectSEDlike, Data=Data, method = 'BFGS')
In practice we might want to use something better like CMA to get pretty close pretty quickly (you can get this from my GitHub at asgr/cmaeshpc). Here we set the wall time to 2 minutes (which is usually enough for this many parameters).
set.seed(1)
library(cmaeshpc)
Data$fit = 'CMA'
badpar = (Data$intervals$lo + Data$intervals$hi) / 2 #CMA is pretty tolerant of terrible initial guesses, unlike optim and LD.
CMAout = cmaeshpc(par=badpar, fn=ProSpectSEDlike, Data=Data, lower=Data$intervals$lo,
upper=Data$intervals$hi, control=list(trace=TRUE, maxwalltime=2))
print(CMAout$par)
There Can Be Only One
Since writing the above example, I have developed the very easy to use Highlander package (available on my GitHub as asgr/Highlander). This combines the above CMA and LD steps painlessly, alternating between CMA and LD phases (by default it does both twice, with 100 iterations for each of the 4 phases) and deals with all the limits auto-magically. The above example can be run simply with:
library(Highlander)
Highout = Highlander(badpar, Data=Data, likefunc=ProSpectSEDlike, Niters=c(1e3,1e3))
This should find a better solution faster than just running CMA or LD on their own. The main arguments you might want to play with would be the number of optimisation phases optim_iters (where a phase is a CMA followed by LD), the number of iteration steps within each process Niters (vector argument specifying CMA followed by LD iterations), and maybe the number of iterations in the very last LD optimisation NfinalMCMC (since you might want to use these chains for posterior exploration).
Exploring Results
We can check the posterior monitored distribution of stellar mass formed (masstot). This is a derived parameter given the star formation history, hence we have to monitor it as we go. Otherwise we would need to reconstruct the mass formed for all 10^4 combinations of parameters generated (which would take a long time!).
maghist(LDout$Monitor[,"masstot"], verbose = FALSE, xlab='Stellar Mass / Msol', ylab='PDF')
abline(v=genSED$Stars$masstot, col='red')
We can also see what fluxes were generated by our various MCMC samples:
magplot(flux_input$pivwave, LDout$Monitor[1,4:23], type='l', log='xy', grid=TRUE,
xlab="Wavelength (Ang)", ylab='Flux Density / Jy')
for(i in 2:1e4){
lines(flux_input$pivwave, LDout$Monitor[i,4:23], col=hsv(alpha=0.1))
}
points(flux_input[,c("pivwave","flux")])
magerr(flux_input$pivwave, flux_input$flux, ylo=flux_input$fluxerr)
ProSpectSED also provides a couple of high level convenience plotting functions:
Data$fit = 'check' #we have to set the fit type to 'check' to get the required outputs
outputs = ProSpectSEDlike(parm=Highout$parm, Data=Data)
plot(outputs)
plot(outputs$SEDout)
Fitting a Galaxy for Redshift
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. In this example we will do something very similar to the above, but now also fit for the redshift. It is easy to add redshift as an explored parameter, given access to SED based photo-z. Basically we have to remove it from arglist and add it to parm.names, making sure we give it a reasonable range. Here we explore it in log-space:
Note since we are fitting for redshift now we can no longer pass LumDist_Mpc because this will change for every realisation.
Data_z=list(flux=flux_input,
arglist=list(massfunc=massfunc_snorm_trunc, agemax=agemax, magemax=agemax/1e9, Z=Zfunc_massmap_lin), # we remove the explicit z from here
speclib=BC03lr,
Dale=Dale_NormTot,
filtout=filtout,
SFH=SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
parm.names=c('z','mSFR','mpeak','mperiod','mskew','tau_birth','tau_screen', # we add z as a parameter to fit here
'alpha_SF_birth','alpha_SF_screen'), # which parameters to fit for
logged=c(T,T,T,T,F,T,T,F,F), # fit parameters in logged or linear space (note z will be fitted in log space)
intervals=list(lo=c(-3,-4,-2,-1,-0.5,-2.5,-2.5,0,0), hi=c(1,3,1,1,1,1.5,1,4,4)), # fitting range for parameters (z will be 0.001 to 10)
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='norm',
verbose=FALSE
)
Now we can run an MCMC fit using LaplacesDemon (depending on machine, this might take 10-20 minutes):
LDout = LaplacesDemon(Model=ProSpectSEDlike, Data=Data_z, Initial.Values=c(0.1,inpar),
control=list(abstol=0.1), Iterations=1e4, Algorithm='CHARM', Thinning=1)
We can see that the redshift, fitted without any prior (or implicitly a uniform in log-space prior) is degenerate with many other parameters. In particular the total stellar mass, where higher redshifts results in more stellar mass. This is because a galaxy need to be brighter and on average more massive if placed at higher redshift to create the same amount of observed flux.
magtri(cbind(LDout$Posterior1,SM=LDout$Monitor[5001:1e4,"masstot"]))
In practice we might want to use something better like CMA to get pretty close pretty quickly. Here we set the wall time to 2 minutes (which is usually enough for this many parameters).
library(cmaeshpc)
Data_z$fit = 'CMA'
badpar = (Data_z$intervals$lo + Data_z$intervals$hi) / 2 #CMA is pretty tolerant of terrible initial guesses, unlike optim and LD.
CMAout = cmaeshpc(par=badpar, fn=ProSpectSEDlike, Data=Data_z, lower=Data_z$intervals$lo,
upper=Data_z$intervals$hi, control=list(trace=TRUE, maxwalltime=2))
print(CMAout$par)
With photo-z like prior
prior = function(parm){
return(dnorm(parm[1], mean=-1, sd=0.02, log=TRUE)) # 10^0.02 gives about 4.7% photo-z error
}
Data_z2=list(flux=flux_input,
arglist=list(massfunc=massfunc_snorm_trunc, agemax=agemax, magemax=agemax/1e9, Z=Zfunc_massmap_lin), # we remove the explicit z from here
speclib=BC03lr,
Dale=Dale_NormTot,
filtout=filtout,
SFH=SFHfunc, # the preferred functional form of the SFH (eg either SFHfunc, SFHburst)
parm.names=c('z','mSFR','mpeak','mperiod','mskew','tau_birth','tau_screen', # we add z as a parameter to fit here
'alpha_SF_birth','alpha_SF_screen'), # which parameters to fit for
logged=c(T,T,T,T,F,T,T,F,F), # fit parameters in logged or linear space (note z will be fitted in log space)
prior=prior,
intervals=list(lo=c(-3,-4,-2,-1,-0.5,-2.5,-2.5,0,0), hi=c(1,3,1,1,1,1.5,1,4,4)), # fitting range for parameters (z will be 0.001 to 10)
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='norm',
verbose=FALSE
)
LDout_prior = LaplacesDemon(Model=ProSpectSEDlike, Data=Data_z2, Initial.Values=c(0.1,inpar),
control=list(abstol=0.1), Iterations=1e4, Algorithm='CHARM', Thinning=1)