In my opinion the easiest way to do the fit is to use the find_fit function, which is standard in SAGE. For example:
sage: x = [1, 3, 4]
sage: y = [1.2, 4.5, 3.6]
sage: var('a,b,t')
sage: model(t)=a*t+b
sage: data=zip(x,y)
sage: fit=find_fit(data,model,solution_dict=True) # Dictionary is optional but handy
sage: model.subs(fit) # Show the model
t |--> 0.92142857142840007*t + 0.64285714285636364
sage: plot(model.subs(fit),(t,0,5))+points(data,size=20,color='red') # Nice plot
When you have a huge number of points and you want just a polynomial fit, I found that it is (numerically) better to use the polyfit function from numpy:
sage: import numpy as np
sage: a,b=np.polyfit(x,y,1)
sage: a,b
(0.92142857142857137, 0.64285714285714401)