Skip to content

Instantly share code, notes, and snippets.

Embed
What would you like to do?
Fit to an expontial sum in pymc
# 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)
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment