The first difficulty in trying to fit an exponential distribution to data is that the data may not follow an exponential distribution. Nothing grows exponentially forever. Eventually growth slows down. The simplest way growth can slow down is to follow a logistic curve, but fitting a logistic curve has its own problems, as detailed in the previous post.
Suppose you are convinced that whatever you’re wanting to model follows an exponential curve, at least over the time scale that you’re interested in. This is easier to fit than a logistic curve. If you take the logarithm of the data, you now have a linear regression problem. Linear regression is numerically well-behaved and has been thoroughly explored.
There is a catch, however. When you extrapolate a linear regression, your uncertainly region flares out as you go from observed data to linear predictions based on the observed data. Your uncertainty grows linearly. But remember that we’re not working with the data per se; we’re working with the logarithm of the data. So on the original scale, the uncertainty flares out exponentially.
I would view the scale-free behavior in an exponential regression as a feature rather than a bug. It not only scales proportional to the predicted value as it grows, it scales down as it shrinks.
In a regular linear regression, we have y[n] = alpha + beta * x[n] + epsilon[n], where epsilon[n] ~ normal(0, sigma). Rearranging, that gives us the direct sampling y[n] ~ normal(alpha + beta * x[n], sigma), which is how you typically code the model.
If we want to move to an exponential model, there are two ways to do it. A simple approach is to exponentiate the linear predictor to get the non-linear additive regression y[n] ~ normal(exp(alpha + beta * x[n]), sigma). This immediately runs into the problem that it can produce negative values of y. The noise scale is constant rather than proportional to the value.
So what we usually do is instead do a linear regression z[n] ~ normal(alpha + beta * x[n], sigma), then set y[n] = exp(z[n]). Plugging in z[n], we get exp(z[n]) = exp(alpha + beta * x[n] + epsilon[n]) = exp(alpha + beta * x[n]) * exp(epsilon[n]). In other words, now we have multiplicative noise that scales proportional to the value of exp(alpha + beta * x[n]). We have also solved the problem of never getting negative values for y.