Line fitting for synthetic spectra using pymc3

Questions
1

Hello!

I’m trying to write a package to fit gaussian/voigt profiles to synthetic noisy absorption spectra I’ve generated. Here’s the code I’m using:

#mcmc model

def gauss(dv,b):
sqrtpi = 1. / (np.sqrt(np.pi) * b * 1e5)
G = sqrtpi * np.exp(-(dv/b)**2)
return G

def Absprofile(dv,N,b,f,lambda_0):
sigma_0 = r_e * c * (10**N) * np.pi
tt = sigma_0 * f * lambda_0 * 1e-8
tau = tt * gauss(dv,b)
return np.exp(-tau)

    model = mc.Model()
    with model:
        for ii in range(1,n+1):
            vars_dict[ii] = {}
            BoundedNormal = mc.Bound(mc.Normal,lower=0.0)
            vars_dict[ii]['v0'] = mc.Uniform('v0_'+str(ii),wavelengths[0],wavelengths[-1])
            vars_dict[ii]['b'] = mc.Uniform('b_'+str(ii),10,50)
            #vars_dict[ii]['b'] = BoundedNormal('b_'+str(ii),mu=25,sd=20)
            vars_dict[ii]['N'] = mc.Uniform('N_'+str(ii),12,20)
        sd = 1. / mc.Uniform('sd',0,1)**2
        for jj in range(1,n+1):
            vv0 = vars_dict[jj]['v0']
            nn = vars_dict[jj]['N']
            bb = vars_dict[jj]['b']
            vv = wav - vv0
            computed_profile *= mc.Deterministic('flux_obs'+str(jj),Absprofile(vv,nn,bb,0.048,1302))
        observation = mc.Normal('obs',mu=fluxx,sd=sd,observed=flx)  
    with model: 
        trace = mc.sample(10000,tune=5000,chains=2,discard_tuned_samples=True, nuts_kwargs=dict(target_accept=0.95))
    pro = np.ones_like(wav)
    for ii in range(1,n+1):
        profile *= Absprofile(wav - trace.get_values('v0_'+str(ii)).mean(),trace.get_values('N_'+str(ii)).mean(),trace.get_values('b_'+str(ii)).mean(),0.048,1302)
    total *= profile

I calculate regions where absorption lines are detected and setting the initial number of lines to the number of local minima in that region. Each line is characterized by 3 parameters : {v0, N, b} which I’m trying to fit using pymc. I’m attaching example fitsblend_profile2 that I ran with this code (the blue lines are the original spectra and orange ones are the fits). The fitting isn’t happening as expected and I’m not sure what I’m doing wrong since I’m a pymc3 newbie.

For isolated lines, this model works well for the {v0, N} parmeters and gives b values that are off with a very large standard deviation. As shown in the attached image, this fails quite badly for spectra with multiple/blended/overlapping lines.

Any ideas/suggestions are appreciated and let me know if clarifications are needed.

Thanks!

Edit : I can only attach one image at a time.

2

You can take a look of Fitting a spectra of gaussians which is a similar question.