# What is the best (fastest and accurate) way to obtain PDF from the given dataset?

Consider some list of two coordinates, points, generated by some arbitrary function function1 (which is a pdf up to a normalization):

function[x1_, x2_] =
If[x1 < Pi/(5*2),
Exp[-(x2/20)^2 Cos[5 x1]^4]/(x2^2 + 900.) Sin[x1]*Cos[5*x1]^4, 0];
RandVals =
Join[RandomReal[{10^-5., Pi/10}, {10^7, 1}],
RandomReal[{0.5, 350}, {10^7, 1}], 2];
WeightsComp =
Hold@Compile[{{DataUnweighted, _Real, 2}},
Table[function[DataUnweighted[[i]][[1]],
DataUnweighted[[i]][[2]]], {i, 1, Length[DataUnweighted], 1}],
CompilationTarget -> "C", RuntimeOptions -> "Speed"] /.
DownValues@function // ReleaseHold;
WeightsVals = WeightsComp[RandVals];
points = RandomVariate[edist, 5*10^6];


I want to obtain the reversed operation, i.e. to get function1 (up to normalization) from points if having the latter precomputed. This means that the true function is unknown, but I should restore it in the end.

This is how I do this. First, I logarithmize the data (this is to improve the construction of the smooth distribution):

pointslogcomp = Compile[{{points, _Real, 2}}, Log[points], CompilationTarget -> "C", RuntimeOptions -> "Speed"]
pointslog = pointslogcomp[points]; // AbsoluteTiming
{Logx1min, Logx1max} = MinMax[pointslog[[All, 1]]];
{Logx2min, Logx2max} = MinMax[pointslog[[All, 2]]];


Then, I use SmoothKernelDistribution and PDF:

(*Computing smooth distribution*)
smoothkerneldistrX =
SmoothKernelDistribution[pointslog,
"LeastSquaresCrossValidation", {"Bounded", {{Logx1min,
Logx1max}, {Logx2min, Logx2max}}, "Gaussian"},
MaxMixtureKernels -> All, InterpolationPoints -> 300,
MaxExtraBandwidths -> 0]; // AbsoluteTiming
(*PDF*)
pdf[x1_, x2_] =
If[Exp[Logx1min] < x1 < Exp[Logx1max] &&
Exp[Logx2min] < x2 < Exp[Logx2max],
Evaluate[1/x1 1/x2 PDF[smoothkerneldistrX, {Log[x1], Log[x2]}]],
0];


The method works:

pdf1[x1_] =
If[Exp[Logx1min] < x1 < Exp[Logx1max],
Evaluate[
1/x1 PDF[MarginalDistribution[smoothkerneldistrX, 1], Log[x1]]],
0];
Show[Histogram[points[[All, 1]], 500, "ProbabilityDensity",
ScalingFunctions -> {"Log", "Log"}],
LogLogPlot[pdf1[x1], {x1, 10^-5, Pi/10.}]]


but is slow (the timing for smoothkerneldistrX is ~30 seconds). Is it possible to develop a faster method of obtaining the pdf while keeping the same accuracy?

Edit

This behavior is shown by the pdf obtained with points and not with pointslog (the other code is the same):

• Recommendation: simplify your question as much as possible before posting. Do you really want to include terms like 2*5 instead of 10, or 30^2 instead of 900? Commented May 15, 2023 at 0:33
• A smooth histogram (SmoothHistogram[points[[All, 1]]]) looks fine. Why take the logs of both the horizontal and vertical axes?
– JimB
Commented May 15, 2023 at 1:52
• @JimB : if I use SmoothKernelDistribution to points instead of pointslog, I get a wrong behavior of the resulting pdf at small values of x1. Commented May 15, 2023 at 7:30
• What exactly is the "wrong behavior"? If that is the basis for asking this question, then that should be made explicit.
– JimB
Commented May 15, 2023 at 15:17
• @JimB : please see the figure in the updated question. Commented May 15, 2023 at 15:25

There is no anomaly. You are getting exactly what SmoothKernelDistribution is designed to do.

AbsoluteTiming[skd = SmoothKernelDistribution[Log[points[[All, 1]]]];]
(* {4.87891, Null} *)

smallX1Values = Select[Log[points[[All, 1]]], # < -8 &];
Show[Plot[PDF[skd, x], {x, -11, -8}, PlotRange -> {{-11, -8}, {0, 0.00001}},
Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"x1", "Probability density"},
PlotRangeClipping -> False],
ListPlot[Transpose[{smallX1Values, ConstantArray[0, Length[low]]}], PlotStyle -> Red]]


One sees the contributions of the individual kernels centered over each data point. As there are only 21 values of x1 below -8 one can't expect too much as the bandwidth is chosen based on all of the data.

Sometimes using an "adaptive kernel" might help as it spreads out the kernels more where there is a lower concentration of data points and makes the kernels less spread out where there is a higher concentration of points but that doesn't necessarily work out well for your example.

(* Adaptive kernel with arbitrary parameters *)
s = 0.5;
bw = 0.075;

(* Small values of x1 *)
smallX1Values = Select[Log[points[[All, 1]]], # < -8 &];
Show[Plot[PDF[skdAdaptive, x], {x, -11, -8}, PlotRange -> {{-11, -8}, {0, 0.00001}},
Frame -> True, FrameLabel -> (Style[#, Bold, 18] &) /@ {"x1", "Probability density"},
PlotRangeClipping -> False],
ListPlot[Transpose[{smallX1Values, ConstantArray[0, Length[smallX1Values]]}], PlotStyle -> Red]]


(* All values of x1 *)
Show[Histogram[Log[points[[All, 1]]], "FreedmanDiaconis", "PDF"],
Plot[PDF[skdAdaptive, x], {x, -11, Max[Log[points[[All, 1]]]]}, PlotRange -> All, Frame -> True,
FrameLabel -> (Style[#, Bold, 18] &) /@ {"x1", "Probability density"},
PlotRangeClipping -> False]]


There is no need to obtain estimates of the marginal pdf's from the estimated joint pdf. Just fit the marginals directly with the variable of interest (Log[points[[All,1]]] or points[[All,1]].

This is an extended comment mainly to show my confusion as to what you want and what you've done.

It appears that you want to find the marginal distribution of $$X_1$$ from the joint distribution of $$X_1$$ and $$X_2$$. The joint distribution is proportional to

$$\frac{\sin (x_1) \cos ^4(5 x_1) \exp \left(-\left(\frac{x_2}{20}\right)^2 \cos ^4(5 x_1)\right)}{x_2^2+900}$$

where $$10^{-5} \leq x_1 \leq \pi/10$$ and $$1/2 \leq x_2 \leq 350$$.

You go through an interesting approach to get a large sample from that bivariate distribution to obtain a nonparametric estimate of the joint probability density which is then used to get a marginal probability density for $$X_1$$. But if the objective is to obtain an accurate and speedy approximation of the marginal distribution of $$X_1$$, none of that is necessary.

One can use NIntegrate to obtain an accurate approximation of the marginal distribution.

First, obtain the constant of integration for the formula you give to get a legitimate bivariate probability density function.

c = NIntegrate[Exp[-(x2/20)^2 Cos[5 x1]^4]/(x2^2 + 900) Sin[x1]*Cos[5*x1]^4,
{x1, 10^-5, Pi/10}, {x2, 1/2, 350}]
(* 0.000186902 *)


Now define a function that gives the marginal probability density for $$X_1$$:

pdfx1[x1_] := NIntegrate[Exp[-(x2/20)^2 Cos[5 x1]^4]/(x2^2 + 900) Sin[x1]*Cos[5*x1]^4/c,
{x2, 1/2, 350}]
Plot[pdfx1[x1], {x1, 10^-5, Pi/10}]


If you need to take a look at the pdf at small values of x1, that works fine:

Plot[pdfx1[x1], {x1, 10^-5, 10^-4}]


I know that some folks here have no trouble plotting a log scale for the density but I consider that blasphemy.

• Thanks for the answer! However, in reality, I am dealing with precomputed points and I do not know the true distribution (say, I compute the phase space of decay products of some particle at the laboratory inertial frame of the latter; and I start with the dataset with momenta of the decaying particles). I have to obtain the pdf from points. I apologize for not stating this clearly in the question. Commented May 15, 2023 at 20:28
• Good. If you could clarify that in your question, then I think you'd get better answers. Also, it appears you want a univariate pdf so why not start with a univariate set of points? Getting an estimate of a bivariate distribution only to want a pdf for a marginal is totally unnecessary. Finally, "I want to obtain the reversed operation" is not possible if you mean figuring out the functional form of the underlying joint pdf unless you know the structure of the joint pdf and only a few parameters are unknown.
– JimB
Commented May 15, 2023 at 21:33