I have some experimental data (vectors x,y) and I've estimated the area under the curve via the trapezoid method (A = trapz(y,x))
The question is, how can I determine the error of this calculation? Is there a Matlab function that can estimate it?
Thanks a lot..
No products are associated with this question.
As Matt J has indicated with his reference to
the error can be estimated in terms of the second derivative of your data function and the widths of your trapezoidal intervals. You can estimate the second derivative in terms of the typical second finite differences in the data divided by the square of the interval widths.
I'm going to revive this just because I think this problem is fairly widespread and there isn't much available by way of resources.
When discussing the error of trapz, it is most commonly understood to be referring to the difference between the integral of an actual underlying function and the estimate made using trapz. That is not the issue here. When working with experimental data, there is no known underlying function. The problem is that the data points themselves are unreliable. The analogous case would be, if you had a known function, every time you call it there is a random error term added to it. The question is to obtain the uncertainty in the integrated area given the uncertainty in each of the data points.
There are a couple of approaches one could take. There are some formulas available here http://cmd.inp.nsk.su/old/cmd2/manuals/cernlib/shortwrups/node88.html, http://wwwasdoc.web.cern.ch/wwwasdoc/shortwrupsdir/d108/top.html
My solution was to implement the trapz() algorithm by hand, and to manually take care of the error propagation at each step. That is, you know how to convert a pair of X and Y values into an estimate of the area for a single segment, and you can use error propagation http://en.wikipedia.org/wiki/Propagation_of_uncertainty (check the table) to work out the uncertainty in that segment. You can then continue propagating the errors as you add segments together. Here is my example code (where x and y have 2 columns, and the second column is the standard deviation for each data point).
int = 0; int_err = 0; for j = 1:length(x)-1 yterm = 0.5*(y(j+1,1)+y(j,1)); xterm = (x(j+1,1)-x(j,1)); yerr = sqrt(0.5*(y(j,2)^2+y(j+1,2)^2)); xerr = sqrt(x(j,2)^2+x(j+1,2)^2); z = yterm * xterm; zerr = sqrt(z^2*((yerr/yterm)^2 + (xerr/xterm)^2)); if ~isfinite(zerr) zerr = 0; end int = int + z; int_err = sqrt(int_err^2 + zerr^2); end