# Estimating error when bootstrapping results are skewed

I'm trying to figure out the most robust way to estimate the error on the best-fit parameters of an exponential model ($$y = x^a$$) for some data. I have heard that bootstrapping is a solid way to attempt this and so have fit models to 1000+ random samples of the data (generated with replacement) and saved the best-fit $$a$$ for each sample.

The problem is that the distribution of these bootstrapped values looks extremely skewed and is nowhere close to Gaussian. I know that with sample statistics like mean, standard error, etc the Central Limit Theorem would imply that the distribution of results approximates a Gaussian as the sample size increases. I don't think this would apply to model parameters, so I suppose it is not super surprising that the bootstrapped distribution is not.

I'm not sure how to express the uncertainty on my best-fit estimate of $$a$$ since typical measures like standard error wouldn't apply. Is there some sort of established method to estimate this sort of error on the peak of a skewed distribution? Thank you for any help.

• "the distribution of results approximates a Gaussian as the sample size increases." Note that the 'sample size' that's most relevant is your original data sample size, not the number of bootstrap samples you evaluated. Also, some statistics are inherently biased and skewed; see this page for example.
– EdM
Sep 17, 2020 at 12:57

(1) The Bootstrapped Distribution: the simplest way to express the uncertainty in $$a$$ would be to simply plot the bootstrapped distribution of values of $$a$$ using a histogram, and let your audience draw their own conclusions about the relationship between x and y. This won't allow you to do any hypothesis testing around $$a$$, but if your goal is just to express the precision of your estimate then it's perfectly acceptable.
(3) Transform the data: from the way you've phrased your question it sounds like you've been estimating $$a$$ by with something like $$\hat{a} = \text{argmin}_a \sum_i |y_i - x_i^a|^2$$, where $$i$$ is indexing over your datapoints. This isn't wrong per se, but it is inconvenient for the reasons you're encountering. If you're willing to change the loss function, it would be far simpler to simply fit the model $$\log(y) = a \log(x)$$. As you noted, the coefficient estimates of a linear regression are normally distributed, so by transforming your data you simplify the problem enormously, since now you have simple linear regression (with the minor twist the coefficient is know to be equal to 0, which doesn't matter for the results you need). This is the approach I would recommend, assuming it doesn't violate the constraints of your problem for some other reason.