0
$\begingroup$

Consider a dataset data-angles-gamma-factors.txt having the form {{x1,y1},{x2,y2},...}. Each row represents an event with values x,y; all events have the equal weight. I would like to make a double distribution function using this dataset. I do the following:

SetDirectory[NotebookDirectory[]];
Data1 = Import["data-angles-gamma-factors.txt", "Table"];
DoubleDistributionTemp = 
  SmoothKernelDistribution[Data1, MaxExtraBandwidths -> 0];

To check the correctness of the double distribution, I check the single distribution for x obtained by integrating the double distribution and the single distribution obtained directly using the initial dataset:

BinsNumber = 300;
BinsHeight = IntegerPart[Length[Data1]/BinsNumber];
binsx = 
  Partition[Sort[#], Round[Length[#]/#2]] &[Data1[[All, 1]], 
   BinsNumber];
xDistrTable = 
  Table[{(binsx[[i]][[1]] + binsx[[i]][[BinsHeight]])/
    2, BinsHeight/(
    binsx[[i]][[BinsHeight]] - binsx[[i]][[1]])}, {i, 
    1, BinsNumber - 1, 1}];
xDistr[x_] = 
 Quiet[Interpolation[xDistrTable, 
    InterpolationOrder -> 0][x]/
  NIntegrate[
   Interpolation[xDistrTable, 
     InterpolationOrder -> 0][X], {X, 0, Pi}]]
ListLogPlot[{Table[{x, xDistr[x]}, {x,
     0, Pi, 10^-2}], 
  Table[{x, 
    NIntegrate[
     PDF[DoubleDistributionTemp, {x, y}], {y, 
      1, 41.9}]}, {x, 0, Pi, 10^-2}]}, 
 PlotStyle -> {Blue, Red, Darker@Darker@Green}, Frame -> True, 
 FrameStyle -> Directive[24, Black], 
 PlotRange -> {{0, Pi}, {10^-2, 10}}, ImageSize -> Large, 
 GridLines -> Automatic]

However, the results are different, which tells that the double distribution function is incorrect:

enter image description here

Could you please tell me how to perform the interpolation correctly?

$\endgroup$
1
$\begingroup$

I'm assuming you want an estimate of a univariate density from a bivariate density using SmoothKernelDistribution. Both the x and y values are bounded so with ` you need to use the "Bounded" option and because the data is nothing like a bivariate normal distribution, a better choice for the bandwidth is necessary.

First I'll illustration things with a univariate approach and then move on to the original question about getting an appropriate/acceptable univariate density estimate from a bivariate density estimate.

The first coordinate is bounded by 0 and $\pi$. Below is the code to estimate the univariate density for the first coordinate using three different approaches to setting the bandwidth: Automatic (the default based on a univariate normal), a fixed value (0.002 in this case), and through least squares cross-validation (which in general has the best track record).

Data1 = Import[data-angles-gamma-factors.txt", "Table"];

(* Estimate univariate density with 3 different approaches for setting the binwidth *)
oneDStdGaussian = SmoothKernelDistribution[Data1[[All, 1]], Automatic,
   {"Bounded", {0, π}, "Gaussian"}, MaxMixtureKernels -> All, MaxRecursion -> 6];
oneD002 = SmoothKernelDistribution[Data1[[All, 1]], 0.002,
   {"Bounded", {0, π}, "Gaussian"}, MaxMixtureKernels -> All, MaxRecursion -> 6];
oneDLSCV = SmoothKernelDistribution[Data1[[All, 1]], LeastSquaresCrossValidation",
   {"Bounded", {0, π}, "Gaussian"}, MaxMixtureKernels -> All, MaxRecursion -> 6];

(* Create histograms for the background *)
h = Histogram[Data1[[All, 1]], 1000, "PDF", PlotRange -> {{0, π}, Automatic}, Frame -> True];
hLeft = Histogram[Data1[[All, 1]], 1000, "PDF", PlotRange -> {{0, 0.1}, Automatic},
   Frame -> True, PlotRangeClipping -> True];
hRight = Histogram[Data1[[All, 1]], 1000, "PDF", PlotRange -> {{π - 0.1, π}, Automatic},
   Frame -> True, PlotRangeClipping -> True];

(* Combine graphics *)
Show[h,
 Plot[{PDF[oneDStdGaussian, x], PDF[oneD002, x], 
   PDF[oneDLSCV, x]}, {x, 0, π}, PlotRange -> All,
  PlotLegends -> Placed[{"Automatic", "0.002", "LSCV"}, {Center, Center}]]]
Show[hLeft, 
 Plot[{PDF[oneDStdGaussian, x], PDF[oneD002, x], PDF[oneDLSCV, x]}, {x, 0, 0.1},
  PlotStyle -> Thickness[0.015],
  PlotLegends -> Placed[{"Automatic", "0.002", "LSCV"}, {Scaled[{0.70, 0.6}], {0, 0}}]]]
Show[hRight, 
 Plot[{PDF[oneDStdGaussian, x], PDF[oneD002, x], PDF[oneDLSCV, x]}, {x, π - 0.1, π},
  PlotStyle -> Thickness[0.015],
  PlotLegends -> Placed[{"Automatic", "0.002", "LSCV"}, {Scaled[{0.1, 0.6}], {0, 0}}]]]

Density estimates with histogram

Left-tail of density estimates

Right-tail of density estimates

Looking at the left and right tails of the histogram one sees a bit of avoidance of values near 0 and $\pi$, respectively. The default bandwidth doesn't do well with this data. The least-squares bandwidth (something slightly larger than 0.002) seems OK and maybe a bandwidth of 0.002 seems a little too bumpy. But that's your call.

So, in short, the "Bounded" option is needed and an appropriate setting for the bandwidth.

Estimating a marginal density from a bivariate density

For estimating the univariate density of the first coordinate one also needs the "Bounded" option for both the first and second coordinate. The second coordinate is bounded by 1 and $\infty$. The estimate of the bivariate density can be found with the following:

twoD = SmoothKernelDistribution[Data1, "LeastSquaresCrossValidation",
   {"Bounded", {{0, π}, {1, ∞}}, "Gaussian"}, 
   MaxMixtureKernels -> All, InterpolationPoints -> 1000];

ContourPlot[PDF[twoD, {x, y}], {x, 0, π}, {y, 1, 20}, 
 PlotPoints -> 100, PlotRange -> All]

Bivariate density estimate

Rather than doing brute force integration Mathematica has the MarginalDistribution function which does all of the integration for you. Here's that density along with the histogram from above:

marginal = MarginalDistribution[twoD, 1];
Show[h, Plot[PDF[marginal, x], {x, 0, π}, PlotRange -> All]]

Marginal density function with histogram in the background

This approach does have a much better fit when using the "Bounded" option and the "LeastSquaresCrossValidation" bandwidth selection option. However, if you really need the univariate density estimate, I would stick with starting just with the data from the first coordinate.

$\endgroup$
  • $\begingroup$ Thank you! The method you shown gives better matching. However, this double distribution (and even single distribution) you proposed does not match the data completely. For an example, there is a discrepancy for small (and large) values of the first coordinate. Namely, the pdf gives 0.0039 for the integral over all y and x from 0 to 1/480, while the data (using Select) gives something like 0.0011. Can this be fixed somehow? $\endgroup$ – John Taylor Jun 23 at 16:53
  • 1
    $\begingroup$ That's part of what I'll address in my updated answer. (Too many chores to do first.) Note that there is a tendency for whatever instrument you use that avoids values near 0 and $\pi$ for the first coordinate. Consider the figures from the following: Histogram[Data1[[All, 1]], 1000, "PDF", PlotRange -> {{0, 0.1}, Automatic}, Frame -> True], Histogram[Data1[[All, 1]], 1000, "PDF", PlotRange -> {{\[Pi] - 0.1, \[Pi]}, Automatic}, Frame -> True]. $\endgroup$ – JimB Jun 23 at 17:31
  • $\begingroup$ Thank you again! However, unfortunately, the obtained function still overestimates the domain of small x. Instead of 0.0039 (for the old distribution) it gives 0.0029. I tried to integrate only from the domain for which there is the available data, but the result did not become sufficienctly better. Should I replace the minimal value of x in the options of the SmoothKernelDistribution by the minimal allowed value from the data? $\endgroup$ – John Taylor Jun 24 at 6:23
  • 1
    $\begingroup$ Got it. So...I'm thinking that there might be 3 ways to deal with this: (1) Don't deal with it, (2) Choose small enough bandwidths for the x and y until there is a closer match, and (3) transform both the x and y to distributions that don't have a very influential boundary issue such as ArcSin[Sqrt[Data1[[All, 1]]/\[Pi]]] and Log[Log[Data1[[All, 2]]]], respectively. I'll try (2) and (3) when I get back in about 8 hours. $\endgroup$ – JimB Jun 24 at 15:48
  • 1
    $\begingroup$ I hope I've been of some help. In my (maybe sometimes inappropriate) sermons, I mention that the real objective and restrictions need to be spelled out. From what you've written, maybe it's just the tail regions that are of interest. If so, then just that data should be looked at separately. These methods (SmoothKernelDistribution) go for an "overall" fit rather than concentrating on more specific aspects. Another available option is an "Adaptive" kernel but while theoretically that would be better, I haven't gotten it to work on even moderately-sized datasets. $\endgroup$ – JimB Jun 24 at 20:54

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.