# Finding the integral of a fitted function

I have a function obtained by fitting some data, and I do not have access to the data itself. The fitting parameters of the function have confidence bounds. I need to obtain an expression for the integral of this function with corresponding errors.

I have thought of solving the problem with a Monte Carlo approach. By integrating perturbations of the fitting function by varying the fitting parameters (assuming some probability distribution, e.g. Gaussian), I could then fit this mock data by a new function, and obtain new fitting parameters with new confidence intervals. Is there a more general procedure, or is this a decent approach? This approach assumes no covariance among the fitting parameters which may or may not be a bad assumption. Any good sources on how to approach this problem are greatly appreciated.

EXAMPLE

Suppose we have the fitted function $f(x)=x^{\alpha} \exp\left(- x^{\beta}\right),\alpha=1\pm0.02,\beta=3\pm0.3$. We then wish to obtain an expression for $F(t)=\int^{\infty}_tf(x) dx, t>0$. How do I find the best fit to this function, and how do the errors of the $f$ fit propagate to $F$?

• This is really no different than estimating confidence bounds for any other function of the parameters. The main difference, if any, may be that the correlation among the parameter estimates could play an enormous part in the estimates. Do you at least have access to some information about this correlation? I am also suspicious that this is not really a fit to data, but an estimate of a probability distribution for the data. (The two are very different.) Which is it and how was that estimate obtained? What do you know about the range of data involved? (That, too, is crucial.)
– whuber
Mar 13, 2013 at 23:05
• Incidentally, for $\alpha$ and $\beta$ in these ranges, your integrals diverge, because $x^\alpha$ is nearly linear and $\exp(x^{-\beta}) \ge 1$ for positive $x$. Are you sure the negative signs are correctly placed? (I first read it as $x^\alpha \exp(-x^\beta)$, which does converge.)
– whuber
Mar 13, 2013 at 23:15
• Thank you for your comment. It is indeed a fit to data, and I have no a priori knowledge about correlations among the parameters. The function I chose was just an arbitrary (non trivial) example, I will fix the OP to a more sensible function (and provide constraints on the data). You were right in your observation about $\beta$ - that was a typo. The data was obtained by measurements on a population. Mar 14, 2013 at 0:36

Presume you have, as stated, intervals which contain all possible values of the parameters. You can apply the techniques of interval analysis (interval arithmetic), using for instance, INTLAB, to evaluate all possible values of an integral of that function. Integration in the multivariate case (two or more interval parameters) is more complicated than in the univariate (one interval parameter) case.

An excellent modern introduction to interval analysis is the book Introduction to Interval Analysis, by Moore, Kearfott, Cloud . In particular, chapter 9, "Integration of Interval Functions" addresses what you want to do, but is unfortunately not included in the free Google preview. You can get 30% off buying directly from SIAM if you are a SIAM member. It's a fascinating and powerful field, so you may feel it worth the investment of time and money to learn about it (integration is but one tiny application).

A less easily understood, older treatment of integration of interval functions is available at http://www.nsc.ru/interval/Library/Thematic/General/Rall82.pdf .

I'll discuss it with a fixed $t$. So I'll write $F(\alpha,\beta)=\int^{\infty}_tx^{\alpha} \exp\left(- x^{\beta}\right) dx$.

And I'll assume that "$\alpha=1\pm0.02,\,\beta=3\pm0.3$" means that $\alpha$ and $\beta$ are normal-distributed with means $\mu_\alpha=1$, $\mu_\beta=3$ and standard deviations $\sigma_\alpha=0.02$, $\sigma_\beta=0.3$:

$$\rho(\alpha) = \frac{1}{\sigma_\alpha\sqrt{2\pi}}e^{-\frac{(\alpha-\mu_\alpha)^2}{2\sigma_\alpha^2}},\quad \rho(\beta) = \frac{1}{\sigma_\beta\sqrt{2\pi}}e^{-\frac{(\beta-\mu_\beta)^2}{2\sigma_\beta^2}}$$

From this you can get an expected value for $F$: $$E[F]=\int d\alpha\int d\beta\; \rho(\alpha)\rho(\beta) F(\alpha,\beta)$$ and the second momentum $$E[F^2]=\int d\alpha\int d\beta\; \rho(\alpha)\rho(\beta) \left(F(\alpha,\beta)\right)^2$$ thus getting the variance: $$Var(F) = E[F^2] - E[F]^2$$