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)
```