In asgr/ProSpect: Professional Spectral Analysis Package

Making Mock Data

Load the libraries we will need:

library(ProSpect)
library(LaplacesDemon)
library(foreach)
library(celestial)
library(magicaxis)

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:

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],
                   z=redshift,
                   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)

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=redshift, H0=67.8, OmegaM=0.308)

Data=list(flux=flux_input,
          arglist=list(z=redshift, 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)