find_fit_bug.sage

data = [(0,1),(1,2),(2,3),(3,3),(7,2),(8,5),(9,2),(11,1),(12,2),(13,4),(14,6),(15,11),(16,13),(17,14),(18,13),(19,15),(20,10),(21,3),(22,4),(23,12),(24,17),(25,4),(26,2),(27,5),(28,6),(29,6),(30,5),(31,7),(32,6),(33,2),(34,9),(35,5),(36,1),(37,4),(38,5),(39,1)]
var('a1,b1,c1,a2,b2,c2,x')
model(x) = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2)
sol = find_fit(data,model)
show(sol)
f(x) = model(a1=sol[0].rhs(),a2=sol[1].rhs(),b1=sol[2].rhs(),b2=sol[3].rhs(),c1=sol[4].rhs(),c2=sol[5].rhs())
er = sum([pow(f(p[0])-p[1],2) for p in data])/len(data)
show(er)
p2 = scatter_plot(data)
p3 = plot(f(x),(0,40))
p3 + p2