Modelling problem: y_est = (1-a)*specarr_old + a*np.squeeze(SSParr[np.argmin(abs(SSPage - b))])

Questions
1

I am trying to fit a galaxy spectrum using PYMC3, but it seems to me the following code can not work,
could someone please help me in figuring out what’s wrong in the following code?

with pm.Model() as model:
    a = pm.Uniform('a', 0., 1.) #fburst
    b = pm.Uniform('b', 0., 2.) #burst_age        
    sigma = pm.Uniform('sigma', 0., 2.)
    y_est = (1-a)*specarr_old + a*np.squeeze(SSParr[np.argmin(abs(SSPage - b))])
    y_obs = pm.Normal('y_obs', mu=y_est, sigma=sigma, observed=ydata)
    trace = pm.sample(4000)

In this case, there are two free parameters (a & b), the fitting result is really poor. But if I fix b, and change the code as:

   specarr_young = np.squeeze(SSParr[np.argmin(abs(SSPage - b))])
   y_est = (1-a) * specarr_old + a * specarr_young

then the code works. I am not sure why

   y_est = (1-a)*specarr_old + a*np.squeeze(SSParr[np.argmin(abs(SSPage - b))])

does not work.

By the way, in the code:
SSParr is an array with shape (177, 3276),
SSPage is an array with shape (177,),
ydata is an array with shape (3276,))