Ultimately I would like to do model selection to find the number of atomic emission spectral lines in a measured spectrum. I would first like to be able to specify the number of gaussians to fit to a spectrum. I see there is a normal mixture and maybe thats what I should be using but I can’t make sense of it. I have been able to fit data with two peaks with the but with the code below I would have to make a function every time I want to specify a different number of gaussians. What I want is a way to have a way to specify an arbitrary number of gaussians.
def gaussian(x,amp,cent,sig): return amp*np.exp( -1*(x-cent)**2/2/sig**2) amp = 1 cent = 0 sig = 0.1 x = np.linspace(-2,2,100) f_true = gaussian(x,amp,cent,sig) + gaussian(x,1,-1,sig) noise = np.random.normal(size=len(x))*.002 f = f_true + noise with Model() as gaussian_model: #Parameters with priors amp = Uniform('amplitude',lower=0,upper=2) cent= Uniform('centroid',lower=-2,upper=2) sig = Uniform('sigma',lower=0, upper=4) amp2 = Uniform('amplitude2',lower=0,upper=2) cent2= Uniform('centroid2',lower=-2,upper=2) sig2 = Uniform('sigma2',lower=0, upper=4) sigma = Uniform('sigmat',lower=0,upper=1) gauss = Deterministic('gauss',amp*np.exp( -1*(cent -x)**2/2/sig**2)+amp2*np.exp( -1*(cent2 -x)**2/2/sig2**2)) likeli = Normal('y', mu = gauss, sd=sigma,observed = f) start=find_MAP() step=NUTS() trace=sample(1000,step,start=start)
I tried to do something like the snippet below but this caused a huge slow down as well as a terrible fit.
def multi_gaussian(x,amp,cent,sig): function=0 for i in range(0,len(amp)): function = function + amp[i]*np.exp(-1*(x-cent[i])**2/2/sig[i]**2)