1
$\begingroup$

I have created an empirical distribution that given an array containing the points px (e.g.{px0,px1,px2,px100}) the correspondent probability is given by a liner interpolation of the tipe:

enter image description here

The probability density function (PDF) corresponds to the angular coefficient m of the lines, which is computed in the function below.

GeneralEmpiricalDistribution[xdata_, pt_] := 
 Block[{n, tab, cdf, pdf, sz, i, index, m, x, x0, y0, xf, yf},
  sz = Length[xdata];
  n = sz - 2;
  tab = Table[{xdata[[i]], (i - 1)/(n + 1)}, {i, 1, n + 2}];
  For[i = 1, i < sz, i++,
   If[Between[pt, {xdata[[i]], xdata[[i + 1]]}], index = i; Break[] ];
   If[i == sz - 1, index = sz - 1]
   ];
  {x0, y0} = tab[[index]];
  {xf, yf} = tab[[index + 1]];
  m = (yf - y0)/(xf - x0);
  cdf = m (x - x0) + y0;
  pdf = m;
  {cdf, pdf}
  ]

Then to validate the function above, consider the data (xdata) generated from a normal distribution with mean=2, and sdev=0.4:

mean = 2;
sdev = 0.4;
xdata = Table[Quiet[NSolve[CDF[NormalDistribution[mean, sdev], x] == prob, x][[1,1,2]]], {prob, 0.00000029, 1., 0.00999971}];
n = Length[xdata];
xmax = xdata[[n]];
Plot[{GeneralEmpiricalDistribution[xdata, x][[2]], 
PDF[NormalDistribution[mean, sdev], x]}, {x, 0, xmax}]

enter image description here

Which seems to be OK. But when i try to integrate to find the mean again the solution are very different:

NIntegrate[x GeneralEmpiricalDistribution[xdata, x][[2]], {x, 0, xmax}]
NIntegrate[x PDF[NormalDistribution[mean, sdev], x], {x, 0, xmax}]

Results:

0.0959046

1.99989

Why the integration results are different?

$\endgroup$
  • $\begingroup$ NIntegrate doesn't seem to be holding its first argument, so you're integrating GeneralEmpiricalDistribution[xdata, x]. Which is basically a flat-line. $\endgroup$ – b3m2a1 Sep 11 '17 at 20:08
3
$\begingroup$

I mentioned this in a comment, but NIntegrate is evaluating the first argument, presumably so it can try to speed up the integration. Here's a way around that:

GeneralEmpiricalDistributionPDF[xdata_, pt_?NumericQ] :=
 GeneralEmpiricalDistribution[xdata, pt][[2]]

NIntegrate[x GeneralEmpiricalDistributionPDF[xdata, x], {x, 0, xmax}]
NIntegrate[x PDF[NormalDistribution[mean, sdev], x], {x, 0, xmax}]
NIntegrate::slwcon: Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

NIntegrate::ncvb: NIntegrate failed to converge to prescribed accuracy after 9 recursive bisections in x near {x} = {2.93247}. NIntegrate obtained 1.9980823464466289` and 0.0025316961179251796` for the integral and error estimates.
1.99808

1.99989

Note that it complains about the function converging too slowly (might be because it's step-wise?). But the result is pretty much right.


Speeding up with Interpolation

By the way, here's a way to make the integration much faster and also more precise:

interpolatingGED[xdata_] :=

 With[{range = Range @@ Append[.1]@MinMax@xdata},
  Interpolation@
   Table[{x, GeneralEmpiricalDistribution[xdata, x][[2]]}, {x, 
     range}]
  ]

igedPDF[ig_, x_?NumericQ] :=
 ig[x]

With[{ig = interpolatingGED[xdata]},
 NIntegrate[x igedPDF[ig, x], {x, 0, xmax}] // RepeatedTiming
 ]

{0.016, 1.99893}

Basically we use an InterpolatingFunction to avoid the piece-wise discontinuous issue J.M. points out.


Speeding up with Compile

We can also write a compiled version of this which will be faster than the original.

gedParameters =
  Compile[{
    {xdata, _Real, 1},
    {pt, _Real}
    },
   Block[{n, tab, sz, i, index = -1, m, x0, y0, xf, yf},
    sz = Length[xdata]; n = sz - 2;
    tab = Table[{xdata[[i]], (i - 1)/(n + 1)}, {i, 1, n + 2}];
    Do[
     If[xdata[[i]] <= pt <= xdata[[i + 1]],
      index = i; Break[],
      If[i == sz - 1, index = sz - 1]
      ],
     {i, sz - 1}
     ];
    {x0, y0} = tab[[index]];
    {xf, yf} = tab[[index + 1]];
    m = (yf - y0)/(xf - x0);
    {x0, y0, m}
    ],
    "RuntimeOptions" -> {"EvaluateSymbolically" -> False}
   ];

gedParameters[xdata, 
   RandomReal[{0.0008958563278208521`, 3.6092680307500435`}]] // 
  RepeatedTiming // First

0.000028

GeneralEmpiricalDistribution[xdata, 
   RandomReal[{0.0008958563278208521`, 3.6092680307500435`}]] // 
  RepeatedTiming // First

0.00042

This then speeds up the integration a lot:

Quiet@
  NIntegrate[
   x GeneralEmpiricalDistributionPDF[xdata, x], {x, 0, 
    xmax}] // RepeatedTiming
Quiet@
  NIntegrate[
   x gedParameters[xdata, x][[3]], {x, 0, xmax}] // RepeatedTiming

{0.910, 1.99808}

{0.054, 1.99808}
$\endgroup$
  • $\begingroup$ NIntegrate[] is not very good at piecewise discontinuous functions, unless it has a way to figure where the points of discontinuity are. $\endgroup$ – J. M. will be back soon Sep 11 '17 at 20:17
  • $\begingroup$ @J.M. Good to know. I don't use it much. If the function were made to interpolate, would that make NIntegrate happier? $\endgroup$ – b3m2a1 Sep 11 '17 at 20:20
  • $\begingroup$ Expressing as a Piecewise[] or as an InterpolatingFunction[] ought to suffice, since NIntegrate[] at least knows to split those up in subintervals as needed. $\endgroup$ – J. M. will be back soon Sep 11 '17 at 20:21
  • $\begingroup$ @J.M. I made a compiled version of the OP's function (it's about an order-of-magnitude faster) and gave it a similar holding wrapper for integration but it actually integrates slower. The points here only consider where NIntegrate can intro-spect. It should be able to do so here. Do you have any idea why the un-compiled version is faster? $\endgroup$ – b3m2a1 Sep 11 '17 at 20:25
  • $\begingroup$ That part I'm not sure about, since I don't know what's in NIntegrate[]'s guts. $\endgroup$ – J. M. will be back soon Sep 11 '17 at 20:33
3
$\begingroup$

Since your CDF is piecewise linear, it makes sense to use Interpolation[] on it, with the setting InterpolationOrder -> 1:

xdat = Table[InverseCDF[NormalDistribution[2, 0.4], prob],
             {prob, 0.00000029, 1., 0.00999971}];

cdfEmpirical = Interpolation[MapIndexed[Prepend[(#2 - 1)/(Length[xdat] - 1), #1] &, xdat], 
                             InterpolationOrder -> 1, 
                             "ExtrapolationHandler" -> {Automatic,
                                                        "WarningMessage" -> False}]

where we use the undocumented "ExtrapolationHandler" setting to allow the function to extrapolate without throwing a warning. With this, the PDF is just cdfEmpirical'[x], which we can now use with NIntegrate[]

NIntegrate[x cdfEmpirical'[x], {x, 0, xdat[[-1]]}, 
           Method -> "InterpolationPointsSubdivision"]
   1.99802

where the preprocessor method "InterpolationPointsSubdivision" allows NIntegrate[] to split the function internally into sections NIntegrate[] can easily handle.

$\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.