Currently in scipy and in astropy there are no methods to compute fit parameter errors (apart from some exceptions like leastsq where minimizers do return a covariance matrix).

The simplest method to compute parameter errors is to compute the covariance matrix as the inverse second derivative matrix, which numerically estimates the second derivatives like e.g. here and then computes the inverse matrix.

Question: What is the best API for this?
E.g. pyminuit has a method hesse that estimates the covariance matrix and then exposes the error information in multiple convenient ways to the user:

>>> m.hesse()
>>> m.errors
{'y': 0.99999733060860929, 'x': 2.9999999999998779}
>>> m.covariance
{('y', 'x'): 3.17636386e-16, ('x', 'y'): 3.17636386e-16, ('y', 'y'): 0.99999466122434422, ('x', 'x'): 8.9999999999992681}
>>> m.matrix(correlation=True)
((1.0, 1.0587907808231194e-16), (1.0587907808231194e-16, 1.0))

@nden I know we discussed this, but can you describe the API you have in mind for astropy here?

Side note: actually a function to compute the matrix of second derivatives should be added to scipy.diff and then a convenience function to compute the covariance matrix in scipy.optimize:
http://projects.scipy.org/scipy/ticket/1510
But integrating and maintaining that is quite a bit of work, and we’d have to wait a few years until everyone has it in their scipy, so for now I plan put something like this in astropy.
Nevertheless ... if anyone is interested in contributing this to scipy ... this would be very useful in the long run.