NIntegrate a discrete data set with uncertainties

Question

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?)

Answer 1

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"
    ]
]

Answer 2

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"
    ]
]