[Feature request] uncertainties on best fit parameters

[sefffal]

Documenting this as a feature request, although I might try to implement it myself.

It would seem useful to be able to estimate uncertainties on all PSF model parameters when using fit.

One way this could be implemented today with minimal changes could be to assume Gaussian errors and calculate them using the inverse of the Hessian at the best fitting location. This is a little crude but would be a nice start.

We could also imagine an extension package that uses eg. Turing.jl machinery behind the scenes for Bayesian parameter estimation, but maybe this belongs elsewhere.