Having trouble in spectral fitting
Closed this issue · 6 comments
I am basically trying to fit dwarf galaxies (having mass <10^9 Solar Mass) using MaNGA survey. But not getting a correct fit. Don't know what is happening. No example is also given for spectral fitting. Also the parameter contours are not correct, showing unreasonable results.
import numpy as np
wv1,y1,z1=np.loadtxt("ManGA_12_180451_reduced.txt",unpack=True)
z1 = np.clip(z1, 0.05, np.inf)
from prospect.utils.obsutils import fix_obs
obs = dict(wavelength=wv1,spectrum=y1, unc=z1, redshift=0.278,
maggies=None, maggies_unc=None, filters=None)
obs = fix_obs(obs)
from prospect.models.templates import TemplateLibrary
from prospect.models import SpecModel
model_params = TemplateLibrary["parametric_sfh"]
model_params.update(TemplateLibrary["nebular"])
model_params["zred"]["init"] = obs["redshift"]
model = SpecModel(model_params)
assert len(model.free_params) == 5
print(model)
noise_model = (None, None)
from prospect.sources import CSPSpecBasis
sps = CSPSpecBasis(zcontinuous=1)
print(sps.ssp.libraries)
current_parameters = ",".join([f"{p}={v}" for p, v in zip(model.free_params, model.theta)])
print(current_parameters)
spec, phot, mfrac = model.predict(model.theta, obs=obs, sps=sps)
#print(spec/y1)
from prospect.fitting import lnprobfn, fit_model
fitting_kwargs = dict(nlive_init=400, nested_method="rwalk", nested_posterior_thresh=0.05)
output = fit_model(obs, model, sps, optimize= True, dynesty=True, lnprobfn=lnprobfn, noise=noise_model, **fitting_kwargs)
result, duration = output["sampling"]
from prospect.io import write_results as writer
hfile = "./quickstart_dynesty_mcmc.h5"
writer.write_hdf5(hfile, {}, model, obs, output["sampling"][0], None, sps=sps, tsample=output["sampling"][1],
toptimize=0.0)
from prospect.io import read_results as reader
hfile = "./quickstart_dynesty_mcmc.h5"
out, out_obs, out_model = reader.results_from(hfile)
import matplotlib.pyplot as pl
from prospect.plotting import corner
nsamples, ndim = out["chain"].shape
cfig, axes = pl.subplots(ndim, ndim, figsize=(10,9))
axes = corner.allcorner(out["chain"].T, out["theta_labels"], axes, weights=out["weights"], color="royalblue", show_titles=True)
from prospect.plotting.utils import best_sample
pbest = best_sample(out)
corner.scatter(pbest[:, None], axes, color="firebrick", marker="o")
import matplotlib.pyplot as pl
sfig, saxes = pl.subplots(2, 1, gridspec_kw=dict(height_ratios=[1, 4]), sharex=True)
ax = saxes[1]
pwave = wv1
ax.plot(pwave, out_obs["spectrum"], linestyle="", marker="o", color="k")
#ax.errorbar(pwave, spectr, spec_err, linestyle="", color="k", zorder=10)
bsed = out["bestfit"]
ax.plot(wv1, bsed["spectrum"], color="firebrick", label="MAP sample")
ax = saxes[0]
chi = (out_obs["spectrum"] - bsed["spectrum"]) / out_obs["unc"]
ax.plot(pwave, chi, linestyle="", marker="o", color="k")
ax.axhline(0, color="k", linestyle=":")
#ax.set_ylim(-2, 2)
ax.set_ylabel(r"$\chi_{\rm best}$")
I am doing like this. Is there anything I need to change? I am getting results like this....
Hi, please see some of the discussion at https://prospect.readthedocs.io/en/latest/faq.html#can-i-fit-my-spectrum-too
I suspect you will need at least the "spectral_smoothing" and "optimize_speccal" parameter template sets, e.g. model_params.update(TemplateLibrary["spectral_smoothing"])
I have fitted the spectra successfully. However, the parameters look odd.
Free Parameters: (name: prior)
mass: <class 'prospect.models.priors.LogUniform'>(mini=6309573.44480193,maxi=100000000)
logzsol: <class 'prospect.models.priors.TopHat'>(mini=-0.6,maxi=0.2)
dust2: <class 'prospect.models.priors.TopHat'>(mini=-1,maxi=2.0)
tage: <class 'prospect.models.priors.TopHat'>(mini=1.0,maxi=13)
tau: <class 'prospect.models.priors.LogUniform'>(mini=0.1,maxi=100.0)
sigma_smooth: <class 'prospect.models.priors.TopHat'>(mini=10,maxi=300)
Fixed Parameters: (name: value [, depends_on])
zred: [0.278]
sfh: [4]
imf_type: [1]
dust_type: [0]
lumdist: [125.17384]
imf: [1]
add_dust_emission: [ True]
duste_umin: [1.]
duste_qpah: [4.]
duste_gamma: [0.001]
Initial free parameter vector theta:
[ 1.0e+08 -5.0e-01 5.0e-01 8.7e+00 1.0e+00 2.0e+02]
Initial parameter dictionary:
{'zred': array([0.278]), 'mass': array([1.e+08]), 'logzsol': array([-0.5]), 'dust2': array([0.5]), 'sfh': array([4]), 'tage': array([8.7]), 'imf_type': array([1]), 'dust_type': array([0]), 'tau': array([1]), 'lumdist': array([125.17384]), 'imf': array([1]), 'sigma_smooth': array([200.]), 'add_dust_emission': array([ True]), 'duste_umin': array([1.]), 'duste_qpah': array([4.]), 'duste_gamma': array([0.001])}
The plots:
When including the spectroscopic calibration polynomial you should include at least one photometric popint to set the overall normalization. But also the chain is not even close to converged. You may need to run for many more iterations, and you might also try the dynesty sampling backend.
You may also have to be careful about matching the line spread function of your data, which depends on the instrument.
There is an example for fitting SDSS spectra and photometry as in the prospector paper at https://github.com/bd-j/exspect/blob/main/fitting/psb_params.py
Please reopen if there are further questions, thanks!