0
$\begingroup$

Consider a variable x defined in the range $$ \tag 1 x \in (1,10^{7}) $$

Let us assume that we have three datasets belonging to three different intervals of x:

n1=10^5;
n2=10^6;
n3=10^7;
xmin1=1;
xmax1=10;
xmin2=10;
xmax2=10^4;
xmin3=10^4;
xmax3=10^7;
dataset1=RandomReal[{xmin1,xmax1},n1];
dataset2=RandomReal[{xmin2,xmax2},n2];
dataset3=RandomReal[{xmin3,xmax3},n3];
weights1=Table[(xmax1-xmin1),{i,1,n1,1}];
weights2=Table[(xmax2-xmin2),{i,1,n2,1}];
weights3=Table[(xmax3-xmin3),{i,1,n3,1}];

The weights are chosen such that the combined dataset corresponds to $x$ being uniformly distributed on the whole range $(1)$.

Could you please tell me how, by having the three datasets for three regions, to obtain a smooth distribution function (PDF) for the whole domain?

My attempt is as follows. I restore the distribution per interval in the form

distrtemp[data_, weights_, xmin_, xmax_] := Block[{},
  weighteddat = WeightedData[data, weights];
  edistX = EmpiricalDistribution[weighteddat];
  RandomX = RandomVariate[edistX, 10^7];
  smoothkerneldistr = 
   SmoothKernelDistribution[{Log[#[[1]]]} & /@ Partition[RandomX, 1] //
      Flatten, 
    "LeastSquaresCrossValidation", {"Bounded", {Log[xmin], 
      Log[xmax - xmin]}, "Gaussian"}, MaxMixtureKernels -> All, 
    InterpolationPoints -> 200];
  pdf[x_] = 1/x PDF[smoothkerneldistr, Log[x]];
  {pdf[x], weights[[1]]}]

Having separately the three distributions for three intervals,

distr1[x_] = distrtemp[dataset1, weights1, xmin1, xmax1];
distr2[x_] = distrtemp[dataset2, weights2, xmin2, xmax2];
distr3[x_] = distrtemp[dataset3, weights3, xmin3, xmax3];

an ugly way to "glue" them is to sum with appropriate weights:

distr[x_]=(distr1[x][[2]]*distr1[x][[1]]+distr2[x][[2]]*distr2[x][[1]]+distr3[x][[2]]*distr3[x][[1]])/(distr1[x][[2]]+distr2[x][[2]]+distr3[x][[2]]);

However, this method works badly at the boundaries of each region: instead of a smooth constant line, I get a curve with increasing wiggles:

LogLogPlot[distr[x], {x, 1, 10^7}, PlotRange -> All]

enter image description here

$\endgroup$
1
  • $\begingroup$ I'm going to modify my answer because I think I finally understand what you want. And the fix for your "wiggles" turns out to be extremely simple (and then obvious): You need to have PlotRange -> {Automatic, {0, Automatic}}. In other words the jiggles are relatively small when a zero base is used. $\endgroup$
    – JimB
    Commented Oct 11, 2022 at 20:54

1 Answer 1

3
$\begingroup$

You have a uniform distribution from 1 to 10^7 but you have stratified that range into 3 sections with sample sizes 10^4, 10^5, and 10^7. It is desired to sample the regions separately and then combine everything to estimate the original uniform distribution.

(* Set parameters *)
{n1, n2, n3} = {10^4, 10^5, 10^6}; (* Sample sizes *)
{xmin1, xmax1} = {1, 10}; (* Region boundaries *)
{xmin2, xmax2} = {10, 10^4};
{xmin3, xmax3} = {10^4, 10^7};

(* Generate data *)
dataset1 = RandomReal[{xmin1, xmax1}, n1];
dataset2 = RandomReal[{xmin2, xmax2}, n2];
dataset3 = RandomReal[{xmin3, xmax3}, n3];

(* Estimate distribution for each region *)
 skd1 = SmoothKernelDistribution[dataset1, 
   "LeastSquaresCrossValidation", {"Bounded", {1, 10}, "Gaussian"}];
skd2 = SmoothKernelDistribution[dataset2, 
   "LeastSquaresCrossValidation", {"Bounded", {10, 10^4}, "Gaussian"}];
skd3 = SmoothKernelDistribution[dataset3, 
   "LeastSquaresCrossValidation", {"Bounded", {10^4, 10^7},  "Gaussian"}];

(* Weights are proportional to the range of each region *)
w = {xmax1 - xmin1, xmax2 - xmin2, 
   xmax3 - xmin3}/(xmax1 + xmax2 + xmax3 - xmin1 - xmin2 - xmin3)

(* Overall pdf *)
pdf[x_] := Piecewise[{{w[[1]] PDF[skd1, x], xmin1 <= x < xmax1}, 
  {w[[2]] PDF[skd2, x], xmin2 <= x < xmax2}, 
  {w[[3]] PDF[skd3, x], xmin3 <= x < xmax3}}, 0] // PiecewiseExpand

Plot[pdf[x], {x, -10^6, 1.2*10^7}, PlotRange -> {All, {0, Automatic}},
  Exclusions -> None, PlotRangeClipping -> False]

Overall pdf

Excising the demon wiggles

If one knows that the resulting distribution is uniform, then there would be no need to use nonparametric density estimation (through SmoothKernelDistribution). If one could make statements about the range of possible values of the derivative of the underlying density, then one could appropriately improve on the "look" of the resulting density estimate.

But if only the data is available (along in this case with the bounds), then SmoothKernelDistribution is doing what it is supposed to do. If the criterion for choosing the optimal bandwidth is "how it looks", then the only way to do that is to try different bandwidths.

For these sample sizes, SmoothKernelDistribution won't always perform the requested "LeastSquaresCrossValidation" and falls back to "Silverman". In any event, the wiggles can be toned down by choosing larger bandwidths.

Here is an approach to see the effect of deviating from the bandwidth chosen by SmoothKernelDistribution:

Manipulate[
 skd1x = SmoothKernelDistribution[dataset1, 
   bw1 skd1[[2, 3]], {"Bounded", {1, 10}, "Gaussian"}];
 skd2x = SmoothKernelDistribution[dataset2, 
   bw2 skd2[[2, 3]], {"Bounded", {10, 10^4}, "Gaussian"}];
 skd3x = SmoothKernelDistribution[dataset3, 
   bw3 skd3[[2, 3]], {"Bounded", {10^4, 10^7}, "Gaussian"}];
 
 (*Overall pdf*)
 pdfx[x_] := Piecewise[{{w[[1]] PDF[skd1x, x], xmin1 <= x < xmax1}, 
   {w[[2]] PDF[skd2x, x], xmin2 <= x < xmax2}, 
   {w[[3]] PDF[skd3x, x], xmin3 <= x < xmax3}}, 0] // PiecewiseExpand;
 
 GraphicsColumn[{LogLinearPlot[{pdf[x], pdfx[x]}, {x, 1, 10^7}, 
    PlotRange -> All, Exclusions -> None],
   LogLinearPlot[{pdf[x], pdfx[x]}, {x, 1, 10^7}, 
    PlotRange -> {All, {0, All}}, Exclusions -> None]}],
 {{bw1, 15, "Band width multiplier 1"}, 1, 30, Appearance -> "Labeled"},
 {{bw2, 5, "Band width multiplier 2"}, 1, 30, Appearance -> "Labeled"},
 {{bw3, 12, "Band width multiplier 3"}, 1, 30, Appearance -> "Labeled"},
 TrackedSymbols :> {bw1, bw2, bw3}]

Effect of increasing bandwidth

$\endgroup$
4
  • $\begingroup$ Thank you for the answer! However, the weights are needed to ensure that the total distribution of x in the whole interval from 1 to $10^7$ is uniform. Therefore, the weights for each of the three intervals are the length of the interval. If I would need to merge the three datasets, I would then also divide the weights by the number of entries in datasets. $\endgroup$ Commented Oct 11, 2022 at 6:28
  • $\begingroup$ So, my goal is actually to restore a perfect smooth constant line instead of the wiggly figure I attached in the question. $\endgroup$ Commented Oct 11, 2022 at 7:49
  • $\begingroup$ I've completely changed my answer now that I understand what you wanted to do with the weights. $\endgroup$
    – JimB
    Commented Oct 11, 2022 at 21:12
  • $\begingroup$ Thank you! But in your plot, the region where I had unsmooth behavior cannot be resolved (due tot linear scaling). $\endgroup$ Commented Oct 12, 2022 at 14:17

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.