Fit to an expontial sum in pymc

exp_sum_pymc.py

# Make test data
t=np.linspace(0,20,200)
decays = np.array([2.3,9])
x=np.linspace(-1,1, 32)
amp1 = 10*np.exp(-(x+0.5)**2/0.3)
amp2 = -8*np.exp(-(x-0.5)**2/0.4)
amps=np.column_stack((amp1,amp2))
d=amps.dot(np.exp(-t/decays[:,None]))

#Add Noise
noise1 = np.linspace(0,0.5,32)
np.random.shuffle(noise1)
dn=d+noise1[:,None]*np.random.randn(32,t.size)
dn+=np.random.randn(1, t.size)*0.4

#Define model
sigma=pymc.Uniform('sig',0,1, size=33)
amps=pymc.Uniform('amps',-20,20,size=(32,2))
decays=pymc.Uniform('decays',0,20,size=2)
base =pymc.exp(-t/decays[:,None])
model= pymc.LinearCombination('fit',[amps],[base])
global_err=pymc.Normal('errs',0,1/sigma[-1]**2, size=t.size) 
ch_errs=[]
for i in range(32):
    x=pymc.Normal('err'+str(i), model[i,:]+err_e, tau=1/sigma[i]**2,
                  value=dn[i,:],observed=True,size=t.size)    
    ch_errs.append(x)

                           
@pymc.potential
def orderd(decays=decays):
    """A pair potential"""
    return -10000000*(decays[0]-decays[1]>2)
#
m=pymc.Model([amps, decays, sigma,orderd, model, global_err ]+ch_errs)