# NIntegrate a discrete data set with uncertainties

Consider a data set with some numerical uncertainty in the y-values, such as the following

g[x_] := x^3
data = Table[{x, Around[g[x], RandomReal[{0, .1}]]}, {x, 0, 2, .1}]


I want to numerically integrate the data and find the uncertainty in the result. The only way I know how to numerically integrate the data in Mathematica is to first interpolate, such as

interp = Interpolation[data]
int[y_] := NIntegrate[interp[x], {x, 0, y}]


When I do so with the data with uncertainty, however, I get

The integrand InterpolatingFunction[...] has evaluated to non-numerical values for all sampling points in the region with boundaries {{0,0.1}}

Is Mathematica capable of doing these numerical integrals with uncertainties? And if so, is it documented how it estimates the uncertainty? If not, is there a standard method/workaround to this (I guess something like Monte Carlo resampling?)

• "numerical integrals with uncertainties" - since I am not familiar with how these are carried out in general, out of curiosity could you point to a source of the general method? Feb 14 at 14:50
• Since it doesn't look like Interpolation is making any use of your uncertainties and thus your integral will make no use of your uncertainties, perhaps reference.wolfram.com/language/howto/… might be of interest to you, but I still don't think your uncertainties are going to be directly applicable to NIntegrate
– Bill
Feb 14 at 15:50

# Resampling

Here I ReplaceAll Around with a RandomVariate of the NormalDistribution with the same parameters and make a list with $$1000$$ repeats. For each I do the Interpolation and NIntegrate obtaining a number each time. Then, Around[list] gives an approximate object around the mean of the elements of list and with an uncertainty corresponding to their standard deviation.

int[y_?NumericQ] := Around@Table[
Module[
{
if = Interpolation@ReplaceAll[
data,
Around[a_,b_] :> RandomVariate[ NormalDistribution[a,b] ]
]
},
NIntegrate[ if[x], {x, 0, y}]
]
, {1000}
] Now I compare with the analytical integral $$y^4/4$$

Show[
ListPlot[
Table[
{ k, int[k] }
, {k, 0,2, 0.2}
]
, PlotRange -> {0.0001, All}
, PlotTheme -> "Scientific"
, ScalingFunctions -> "Log"

],
Plot[
y^4/4
, {y, 0,2}
, ScalingFunctions -> "Log"
]
] • Thank you, this is precisely the solution I also came up with
– Kai
Feb 14 at 17:02
• One thing to add to this though, it is worthwhile to check for convergence of the mean and uncertainty as the number of resampled points is varied, although I found 1000 to be more than enough
– Kai
Feb 14 at 17:57

# Using the Trapezoidal_rule

Here the interpolation is happening between the two Nearest points using Rescale[x,{min,max},{ymin,ymax}].

The "Integration" is the trapezium area $$(b-a) \times \left(f(a)+f(b)\right)/2$$.

ClearAll[trap]; (* Trapezoidal integration *)
trap[{a_,fa_},{b_, fb_}]:=(b-a)*(fa+fb)/2

int2[y_]:= Module[
{i,j}, (* Index of the two nearest points *)
{i,j} = Sort@Nearest[ data[[All,1]]->"Index", y, 2];
Rescale[
y,
data[[{i,j},1]],
{#,#+trap@@data[[i;;j]]}&@Total[trap@@@Partition[data[[;;i]], 2,1]]
]
]


Now I compare with the analytical integral.

Show[
ListPlot[
Table[
{ k, int2[k] } (* The function of interest *)
, {k, 0,2, 0.2}
]
, PlotRange -> {0.0001, All}
, PlotTheme -> "Scientific"
, ScalingFunctions -> "Log"
],
Plot[
y^4/4  (* The analytical solution  *)
, {y, 0,2}
, ScalingFunctions -> "Log"
]
] 