Most efficient way to integrate a linear piecewise function

Question

I am having a lot of trouble integrating a simple function. Here's the context: I generate data from two distributions, and I want to calculate the overlap coefficient of representing those datasets with different kinds of histograms. In most cases the resulting PDF is a piecewise linear function, such as a histogram distribution.

These piecewise linear functions are the easiest thing to integrate (high school students could do it by hand), but Mathematica has a real hard time. It's very irregular too. Sometimes is finishes in seconds, sometimes it takes several minutes, frequently it crashes the Kernel. The last one is the real problem because it means I can't collect the data I need.

I've included working code for the basic histogram, but there are many types (16 in total) although some (like the built-in KernelDistribution) are not linear, they are all piecewise. It has crashed at each one of them (some more than others), and it is always on the integration step of each module (yes, in the real code they are Modules). Sometimes the whole thing barely spikes the RAM, other times it uses 32GB, but RAM usage doesn't seem to be correlated with Kernel failure.

The Mathematica documentation says that it automatically considers the boundaries of piecewise for discontinuities, but my hope is that there is some option to set to get it to efficiently calculate these integrals without crashing every few times.

DataPoints1=Sort[RandomVariate[NormalDistribution[RandomReal[{3,7}],RandomReal[{0.5,2}]],100]];
DataPoints2=Sort[RandomVariate[NormalDistribution[RandomReal[{3,7}],RandomReal[{0.5,2}]],100]];
MinDataPoint=Min[Sort[Flatten[{DataPoints1,DataPoints2}]]];
MaxDataPoint=Max[Sort[Flatten[{DataPoints1,DataPoints2}]]];
NumberOfBins=2.*Ceiling[Sqrt[200]];
HistogramBinWidth=(0.002+MaxDataPoint-MinDataPoint)/NumberOfBins;
HistogramBinnedAllData=HistogramList[Sort[Flatten[{DataPoints1,DataPoints2}]],{MinDataPoint-0.001,MaxDataPoint+0.001,HistogramBinWidth}];
HistogramBinBoundaries=HistogramBinnedAllData[[1]];
HistogramBinnedData1=HistogramList[DataPoints1,{HistogramBinBoundaries}][[2]]/(100*HistogramBinWidth);
HistogramBinnedData2=HistogramList[DataPoints2,{HistogramBinBoundaries}][[2]]/(100*HistogramBinWidth);
HistogramPDF1=Piecewise[Join[Table[{HistogramBinnedData1[[i]],HistogramBinBoundaries[[i]]<=x<HistogramBinBoundaries[[i+1]]},{i,1,Length[HistogramBinBoundaries]-2}],{{HistogramBinnedData1[[-1]],HistogramBinBoundaries[[-2]]<=x<=HistogramBinBoundaries[[-1]]}}]];
HistogramPDF2=Piecewise[Join[Table[{HistogramBinnedData2[[i]],HistogramBinBoundaries[[i]]<=x<HistogramBinBoundaries[[i+1]]},{i,1,Length[HistogramBinBoundaries]-2}],{{HistogramBinnedData2[[-1]],HistogramBinBoundaries[[-2]]<=x<=HistogramBinBoundaries[[-1]]}}]];
HistogramPDFOverlap=PiecewiseExpand[Min[HistogramPDF1,HistogramPDF2]];
HistogramDistributionOverlapArea=N@Integrate[HistogramPDFOverlap,{x,MinDataPoint-10,MaxDataPoint+10},Assumptions->x\[Element]Reals]

Here's another complication, that is also maybe a hint to the cause of the problem. When I set NumberOfBins=Ceiling[Sqrt[200]]; the integral usually completes extremely quickly. But when I multiply by two (as in the above code) the time increases by a factor of more than a hundred for the same data. There are at most two times as many integration calculations to perform, so this indicates that the integral function is not doing this in the best way.

Further note: Integrate or NIntegrate is fine with me if one is better, and it doesn't need to be super accurate. I've tried many options and methods based on other posts I've read (TrapezoidalRule seemed good for this, but was much worse), but nothing keeps it from crashing or manages the run time as expected. Any advice or solutions to do this integral properly?

Answer 1

If we know the bin boundaries ahead of time as in the example, the way I would do it is to feed these to NIntegrate[] and skip the PiecewiseExpand simplification. Because of the nature of the pieces of the integrand, a low-order integration rule may be used to save a little speed, but it's optional.

NIntegrate[Min[HistogramPDF1, HistogramPDF2],
  Evaluate[{x, Sequence @@ HistogramBinBoundaries}],
   Method -> {"GaussKronrodRule", "SymbolicProcessing" -> 0, "Points" -> 2}
  ] // AbsoluteTiming
(*  {0.028287, 0.48}  *)

If we bump the number of bins up to NumberOfBins = 10*Ceiling[Sqrt[200]], the timing increases to 0.675531.

Update: It's worth remarking that this approach works best with an "open" rule that avoids the discontinuities in the piecewise function at the bin boundaries. So not the trapezoidal rule, for example, which samples at the end point (bin boundaries). But the midpoint rule would work, if the pdfs are all piecewise constant:

With[{pwints = Partition[
     Union[Cases[HistogramPDF1[[1, All, 2]], _Real, Infinity],
           Cases[HistogramPDF2[[1, All, 2]], _Real, Infinity]],
     2, 1]},
  Dot[
   (Min[HistogramPDF1, HistogramPDF2] /. x -> Mean@# & /@ pwints),
   Flatten[Differences /@ pwints]
   ]
  ] // AbsoluteTiming
(*  {0.002715, 0.48}  *)

With NumberOfBins = 10*Ceiling[Sqrt[200]], the timing increases to 0.061215. So about ten times faster.

Answer 2

Calculating the coefficient of overlapping from samples

Since you are trying to calculate the coefficient of overlapping (OVL) I had a look at what Google has to say and found this paper by Schmid/Schmidt, where nonparametric approaches are described (which I believe is what you are trying to do).

From what I can tell from a cursory look at the paper, some difficulties at finding the integral may be avoided by using a SmoothKernelDistribution. So in your case I - while I am getting some warning - I find this approach to at least get me into the numerical ball park:

RandomSeed[ "09.04.2017" ];

dataPoints1 = RandomVariate[
   NormalDistribution[
       RandomReal[ {3, 7} ], RandomReal[ {0.5, 2} ]
   ], 
   100000 (* using some more data points for clarity *)
];

dataPoints2 = RandomVariate[
   NormalDistribution[
       RandomReal[ {3, 7} ], RandomReal[ {0.5, 2} ]
   ], 
   100000
];

dist1 = SmoothKernelDistribution[ dataPoints1 ];
dist2 = SmoothKernelDistribution[ dataPoints2 ];

NIntegrate[ 
    Min[ PDF[ dist1, x ], PDF[ dist2, x ] ], 
   { x, -∞, ∞ }, 
   WorkingPrecision -> 10
]

0.6315245101

Show @ {
    Plot[ 
        Evaluate[ {PDF[dist1, x], PDF[dist2, x] } ], 
        {x, 0, 10}, 
        Axes -> {True, False}
    ],
    Plot[
        Evaluate[ Min[ PDF[dist1, x], PDF[dist2, x] ] ], 
        {x, 0, 10}, 
        Filling -> Axis
    ]
}

Using HistogramDistribution

From what the OP describes there seem to be issues with regard to using histrogram based PDF. Using HistogramDistribution I cannot detect the described difficulties:

dist1 = HistogramDistribution[ dataPoints1, 2 Ceiling[ Sqrt[200] ] ];
dist2 = HistogramDistribution[ dataPoints2, 2 Ceiling[ Sqrt[200] ] ];

RepeatedTiming @  Quiet @ NIntegrate[
    Min[ PDF[dist1, x], PDF[dist2, x] ], 
    { x, -∞, ∞ }, 
    WorkingPrecision -> 10
]

{5.0297, 0.8113769412}

Using half as many bins will lead to a decrease by a factor of around 25:

dist1 = HistogramDistribution[ dataPoints1, 1 Ceiling[ Sqrt[200] ] ];
dist2 = HistogramDistribution[ dataPoints2, 1 Ceiling[ Sqrt[200] ] ];

RepeatedTiming @ Quiet @ NIntegrate[
    Min[ PDF[dist1, x], PDF[dist2, x] ], 
    { x, -∞, ∞ }, 
    WorkingPrecision -> 10
]

{0.199, 0.8141400000}

Answer 3

I have received many pieces of helpful advice, but not a real answer because combining these things into one method still produces many warnings/errors and sometimes crashes the kernel.

Here is AN answer via updated version of the basic histogram comparison incorporating advice from above. Remember that it is not necessary that the two histograms share the same global bins, although that should make it easier to compute.

First the preliminaries:

DataPoints1=Sort[RandomVariate[NormalDistribution[RandomReal[{3,7}],RandomReal[{0.5,2}]],100]];
DataPoints2=Sort[RandomVariate[NormalDistribution[RandomReal[{3,7}],RandomReal[{0.5,2}]],100]];
MinDataPoint=Min[Flatten[{DataPoints1,DataPoints2}]];
MaxDataPoint=Max[Flatten[{DataPoints1,DataPoints2}]];
NumberOfBins=5*Ceiling[Sqrt[SampleSize1+SampleSize2]];
$MaxPiecewiseCases=(SampleSize1+SampleSize2)*10;

Now the meat of the calculations:

HistogramBinWidth=(0.002+MaxDataPoint-MinDataPoint)/NumberOfBins;
HistogramBinnedAllData=HistogramList[Sort[Flatten[{DataPoints1,DataPoints2}]],{MinDataPoint-0.001,MaxDataPoint+0.001,HistogramBinWidth}];
HistogramBinBoundaries=HistogramBinnedAllData[[1]];
HistogramBinnedData1=HistogramList[DataPoints1,{HistogramBinBoundaries}][[2]]/(SampleSize1*HistogramBinWidth);
HistogramBinnedData2=HistogramList[DataPoints2,{HistogramBinBoundaries}][[2]]/(SampleSize2*HistogramBinWidth);

Using Simplify@ here tends to reduce the number of times that the next step crashes the kernel.

HistogramPDF1=Simplify@Piecewise[Join[Table[{HistogramBinnedData1[[i]],HistogramBinBoundaries[[i]]<=x<HistogramBinBoundaries[[i+1]]},{i,1,Length[HistogramBinBoundaries]-2}],{{HistogramBinnedData1[[-1]],HistogramBinBoundaries[[-2]]<=x<=HistogramBinBoundaries[[-1]]}}]];
HistogramPDF2=Simplify@Piecewise[Join[Table[{HistogramBinnedData2[[i]],HistogramBinBoundaries[[i]]<=x<HistogramBinBoundaries[[i+1]]},{i,1,Length[HistogramBinBoundaries]-2}],{{HistogramBinnedData2[[-1]],HistogramBinBoundaries[[-2]]<=x<=HistogramBinBoundaries[[-1]]}}]];

Here is the biggest problem. In most cases this Simplify@PiecewiseExpand is fine and fast, but it sometimes times out after 300 seconds and crashes the kernel.

HistogramPDFOverlap=Simplify@PiecewiseExpand[Min[HistogramPDF1,HistogramPDF2]];

In some cases the distributions do not overlap, in order to eliminate/reduce the errors of 0 integrals I test for whether the function is zero everywhere explicitly:

HistogramOverlapBoundaries=If[Length[HistogramPDFOverlap]==0,{0.},Union@Cases[HistogramPDFOverlap[[1,All,2]],_Real,Infinity]];
HistogramOverlapCaseValues=If[Length[HistogramPDFOverlap]==0,{0.},Union@Cases[HistogramPDFOverlap[[1,All,1]],_Real,Infinity]];

If the function is not zero everywhere, then do the numerical integration using all the boundaries as exclusions.

HistogramDistributionOverlapArea=If[HistogramOverlapCaseValues=={0.},0,NIntegrate[HistogramPDFOverlap,{x,MinDataPoint-10,MinDataPoint+10},Exclusions->HistogramOverlapBoundaries,AccuracyGoal->8,PrecisionGoal->8,WorkingPrecision->100,MaxRecursion->20]];

When running this code I still get many warnings of the types:

Numerical integration converging too slowly; suspect one of the following: singularity, value of the integration is 0, highly oscillatory integrand, or WorkingPrecision too small.

and a few of the type:

NIntegrate failed to converge to prescribed accuracy after 20 recursive bisections in x near {x} = {3...

Because these are piecewise linear (in fact piecewise uniform) functions, I don't expect to receive errors like this if Mathematica were integrating them appropriately. For example, if it were using the trapezoidal rule with corners at the boundaries, then it should be faster and error free. But I don't know how to get it to it that way. My attempts to get it to use the trapezoidal rule following the documentation made things worse.

So this is AN answer that is better than the version in the question thanks to all your help. It crashes less and runs much faster. However it is not THE answer because I still believe that Mathematica could do this much better given the appropriate conditions.