3
$\begingroup$

I have a list of points which looks like this when plotted with Joined -> True:

enter image description here

I'm interested in the area under the curve with x-axis > 0 as a ratio to the area under the curve as a whole. The obvious way to do this is to define function = Interpolation[list], then use Integrate[function[x],{x,0,Infinity}/Integrate[function[x]],{x,-Infinity,Infinity}].

However, for some of the parameters used to generate this plot, this ratio goes above one. My reading is that the interpolated function goes negative at some point. I know for physical reasons that the function is strictly positive (it is an unnormalized probability). Is there a way to feed this restriction to Interpolation, or perhaps to simply zero out the interpolated function past a certain point?

Edit: some explicit code that shows the problem I'm facing. For this list:

{{-2.,-1.5,-1.,-0.5,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.},{2.447482917*10^-26,7.166525422*10^-23,6.828228711*10^-20,2.221078019*10^-17,2.603521623*10^-15,1.167198034*10^-13,2.134212354*10^-12,1.702246816*10^-11,6.341646627*10^-11,1.183984093*10^-10,1.195443416*10^-10,7.082172125*10^-11,2.655421922*10^-11,6.681130557*10^-12,1.175390799*10^-12}}

(first list is x-axis values, second list is y-axis values). Define interpolating function:

Testfunc = Interpolation[{2.447482916954607`*^-26, 7.166525421661271`*^-23, 6.82822871054717`*^-20, 2.2210780189277698`*^-17, 2.6035216228330743`*^-15, 1.1671980340015243`*^-13, 2.134212354193162`*^-12, 1.7022468159138925`*^-11, 6.341646627292107`*^-11, 1.1839840933620157`*^-10, 1.195443415545412`*^-10, 7.082172125430078`*^-11, 2.655421922388246`*^-11, 6.681130557379854`*^-12, 1.1753907990705456`*^-12}]

Take the numerical integral: NIntegrate[Testfunc[x], {x, 5, 15}]/NIntegrate[Testfunc[x], {x, 1, 15}]

Result is: 1.000008312

|improve this question
$\endgroup$
  • $\begingroup$ Can you provide the data? $\endgroup$ – user21 Sep 9 '19 at 8:20
  • $\begingroup$ @user21 here's one of the lists that show the issue: {{-2.,-1.5,-1.,-0.5,0.,0.5,1.,1.5,2.,2.5,3.,3.5,4.,4.5,5.},{2.447482917*10^-26,7.166525422*10^-23,6.828228711*10^-20,2.221078019*10^-17,2.603521623*10^-15,1.167198034*10^-13,2.134212354*10^-12,1.702246816*10^-11,6.341646627*10^-11,1.183984093*10^-10,1.195443416*10^-10,7.082172125*10^-11,2.655421922*10^-11,6.681130557*10^-12,1.175390799*10^-12}} First list is the x-axis values, second list is the y-axis values. $\endgroup$ – Allure Sep 9 '19 at 8:24
  • $\begingroup$ Can edit your post to contain the data? $\endgroup$ – user21 Sep 9 '19 at 8:25
  • $\begingroup$ Without a better idea of what the function should look like I doubt that you'll be able to get any accurate results. If you look at the curve, the data need to extrapolated significantly to do the integration. And without model, this is not really doable $\endgroup$ – Lukas Lang Sep 9 '19 at 8:44
  • $\begingroup$ @LukasLang I know the function keeps dropping at the extreme ends of the distribution (again for physical reasons), and it shouldn't need to be extrapolated significantly since the values at the ends are very small (note the log scale on the y-axis). That's why when I was doing the numerical integral I only did the integration from x = -2 to 5. $\endgroup$ – Allure Sep 10 '19 at 0:36
6
$\begingroup$

Since the second coordinates are positive you can apply Log and then apply Exp to the function approximating the data before integrating it numerically:

data = {{-2., -1.5, -1., -0.5, 0., 0.5, 1., 1.5, 2., 2.5, 3., 3.5, 4., 4.5, 5.}, 
        {2.447482917*10^-26, 7.166525422*10^-23, 6.828228711*10^-20, 2.221078019*10^-17,
         2.603521623*10^-15, 1.167198034*10^-13, 2.134212354*10^-12, 1.702246816*10^-11,
         6.341646627*10^-11, 1.183984093*10^-10, 1.195443416*10^-10, 7.082172125*10^-11,
         2.655421922*10^-11, 6.681130557*10^-12, 1.175390799*10^-12}};
dataL = Transpose[{data[[1]], data[[2]] // Log}];

f[x_] = Exp[Fit[dataL, x^Range[0,4], x]];

#/(# + #2) &[
  NIntegrate[f[x],{x, 0, Max[dataL[[All,1]]]}],
  NIntegrate[f[x],{x, Min[dataL[[All,1]]], 0}]]

0.999999

Plotting the non-transformed data looks like a normal distribution density, so you could also try:

normal=(E^(-((-m + x)^2/(2 s^2))) n)/(Sqrt[2 Pi] s);
f[x_] = normal /. FindFit[Transpose[data], normal, {m, n, s}, x];

#/(# + #2) &[
  NIntegrate[f[x],{x, 0, Max[dataL[[All,1]]]}],
  NIntegrate[f[x],{x, Min[dataL[[All,1]]], 0}]]

0.999986

|improve this answer
$\endgroup$
3
$\begingroup$

No, you can not specify a limit on the interpolating function. Two options that work for this data set are the use of:

t2=Interpolation[Transpose[data],InterpolationOrder\[Rule]1];

or

t2 = Interpolation[Transpose[data], Method -> "Spline"];
Plot[t2[x], {x, -2, 5}]
NIntegrate[t2[x], {x, 0, 5}]/NIntegrate[t2[x], {x, -2, 5}]
0.9999993190296589`
|improve this answer
$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.