2 edited tags
| link
1
source | link

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