# Most efficient way to integrate a linear piecewise function

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];
HistogramBinWidth=(0.002+MaxDataPoint-MinDataPoint)/NumberOfBins;
HistogramBinnedAllData=HistogramList[Sort[Flatten[{DataPoints1,DataPoints2}]],{MinDataPoint-0.001,MaxDataPoint+0.001,HistogramBinWidth}];
HistogramBinBoundaries=HistogramBinnedAllData[];
HistogramBinnedData1=HistogramList[DataPoints1,{HistogramBinBoundaries}][]/(100*HistogramBinWidth);
HistogramBinnedData2=HistogramList[DataPoints2,{HistogramBinBoundaries}][]/(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]; 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?

• With the same binning for each histogram you can skip all the piecewise functions and just do the sum Total[Min/@Transpose@{HistogramBinnedData1,HistogramBinnedData2}]*HistogramBinWidth. This gives me the same numerical answer as your code, so unless I'm missing something this should work. – N.J.Evans Apr 4 '17 at 13:08
• That works in this case because both histograms use the same global binning from all the data, but that isn't true in general for my different functions. In some cases the min of the two piecewise PDFs will split multiple bars. – Aaron Bramson Apr 4 '17 at 13:14
• Simplify@PiecewiseExpand is very effective here, with that your integral is almost instant. – george2079 Apr 4 '17 at 13:45
• What I would do in that case is create a 0th order interpolation of both binned data sets and create a very fine grained sampling and then use the same method outlined above, assuming george2079's simple fix fails for some case. – N.J.Evans Apr 4 '17 at 13:46
• With v10.1 or later, your definitions for MinDataPoint and MaxDataPoint can be simplified and made more efficient by using MinMax: {MinDataPoint, MaxDataPoint} = MinMax[Sort[Flatten[{DataPoints1, DataPoints2}]]]; Also, it is best to avoid starting user-defined symbols with a capital letter to avoid potential name conflicts with existing or future built-in symbols. – Bob Hanlon Apr 4 '17 at 15:29

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], 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], the timing increases to 0.061215. So about ten times faster.

• That seems to be what the OP is really asking about (+1). :) – gwr Apr 9 '17 at 15:19
• Yes, this is in the right direction. I have found that Nintegrate completes quickly when I feed in the bin boundaries via Exclusions. Although for the histogram the bin boundaries are known, that is not the case in general. In many cases, using Simplify@PiecewiseExpand puts the min function in a format from which I can extract the boundaries. But in some cases Simplify and PiecewiseExpand fail. – Aaron Bramson Apr 9 '17 at 15:19
• @AaronBramson If you have a piecewise function similar to the ones above, something like Union@Cases[First@InternalFromPiecewise[HistogramPDF1], _Real, Infinity] will extract the bin boundaries from each pdf. (If the pdfs are no piecewise-constant and they cross in the middle, NIntegrate might get bogged down and you might need a higher order integration rule -- I'm not sure how well the example reflects all use-cases.) – Michael E2 Apr 9 '17 at 15:31
• @AaronBramson Perhaps simpler: Union@Cases[HistogramPDF1[[1, All, 2]], _Real, Infinity] – Michael E2 Apr 9 '17 at 16:26
• Yeah, that's great. Much simpler than what I had to do the same thing. Getting these boundaries and adding them to exclusions did most of the work. The only thing remaining is a way to detect whether the min function is zero everywhere, and if so to bypass the integral calculation. In a few situations the function is just 0. with only one case, but more often it doesn't simplify down to that. I think Union@Cases[HistogramPDF1[[1, All, 1]], _Real, Infinity] should be the values across all cases...is that right? – Aaron Bramson Apr 10 '17 at 11:57

## 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 ] ];
dist2 = HistogramDistribution[ dataPoints2, 2 Ceiling[ Sqrt ] ];

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 ] ];
dist2 = HistogramDistribution[ dataPoints2, 1 Ceiling[ Sqrt ] ];

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


{0.199, 0.8141400000}

• One can of course play around from here with all the parameters given for SmoothKernelDistribution and NIntegrate. So this answer in this regard is just indicative. – gwr Apr 9 '17 at 14:07
• I am implementing and comparing 16 different nonparametric ways to capture the data (including the kernel distribution and histogram distribution). But most of the techniques are not built into Mathematica, and the built-in versions of others have problems/limitations. The question here is not how to calculate the overlap coefficient based on some data, but rather specifically how to efficiently integrate linear piecewise functions like these. The histogram above is just an easy sample for replication purposes. This is an answer to a different question. – Aaron Bramson Apr 9 '17 at 15:10
• S&S are testing different versions of overlap coefficient estimation using only the Kernel estimate of data. I'm using only the PDF overlap integral measure of the overlap coefficient with data represented in many different nonparametric ways (including Kernel estimation). But they are certainly closely related projects. So I need an integration method that works for all 16 data representations I am using, including custom methods of my own design that are not built into Mathematica. – Aaron Bramson Apr 9 '17 at 15:56
• Minor nit: Show accepts lists, too, so Show@@{...} is unnecessary. – rcollyer Apr 10 '17 at 13:20
• You have to dig. It's in the details section, and oddly in only the last example in the docs ... But, I've seen it used in various places internally, so ... As I said, minor nit. :) – rcollyer Apr 10 '17 at 14:44

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[]; HistogramBinnedData1=HistogramList[DataPoints1,{HistogramBinBoundaries}][]/(SampleSize1*HistogramBinWidth); HistogramBinnedData2=HistogramList[DataPoints2,{HistogramBinBoundaries}][]/(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. • The trapezoidal rule is a poor choice for the reason I explained in my answer. The number of recursive bisections needed to achieve a precision goal$pg$should be at least$\log_2 2\times10^{pg} \Delta y/y$, where$\Delta y$is the difference at the adjacent bin boundaries of a bin whose pdf height is$y$. – Michael E2 Apr 13 '17 at 11:19 • I don't understand what problem you're seeing. On all the examples you've given, the midpoint rule at the end of my post runs much, much faster and gives the same answer as NIntegrate[], if you bump MaxRecursion up above 25, which is the estimate from my previous comment, given your PrecisionGoal and assuming$\Delta y$and$y\$ are between 0 and 1. If something doesn't work, please give an example, and I'll try to address it. (BTW, did you want the upper limit of your NIntegrate to be MinDataPoint+10 or MaxDataPoint+10`?) – Michael E2 Apr 13 '17 at 21:32