We should return confidence intervals on our parameter fits (probably as a model.confidence_ attribute). The standard deviation can be calculated by taking the square root of the diagonal elements in the covariance matrix that is returned by leastsq (np.sqrt(np.diag(pcov)).
It would also be awesome if the predict method could predict the upper and lower confidence interval ranges too!