6
$\begingroup$

I have a sample from distribution whose density is known to be monotonically decreasing. What's the easiest way to incorporate this information into Mathematica for an empirical density plot?

Looking for something like the estimator here which enforces monotonicity by getting density from slopes of the envelope of the CDF.

enter image description here

Example data

ClearAll["Global`*"];
spectralNormalize[vals_] := vals/Max[vals];
empiricalCDF[
   coefs_] := {Most@Sort@coefs, 
    Range[Length@coefs - 1]/Length[coefs]}\[Transpose];
empiricalPDF[
   coefs_] := {Most@Sort@coefs, 
    Most[coefs]*(-1/Differences[coefs])}\[Transpose];

n = 400;
sample := RandomVariate[NormalDistribution[], {n, n}]/Sqrt[n];
vals = SingularValueList[sample];

empiricalCDF = 
 ListPlot[empiricalCDF@vals, PlotLabel -> "empirical CDF"]
empiricalPDF = 
 ListPlot[empiricalPDF@vals, PlotLabel -> "empirical PDF"]

enter image description here

$\endgroup$
3
  • 3
    $\begingroup$ I think this is an excellent (and interesting) question about how to incorporate such "additional information". However, I don't see it as a Mathematica question in that if you knew of an appropriate procedure you wouldn't have any trouble implementing it in Mathematica. I think an answer about a method at either CrossValidated or the Math StackExchange would be where you need to go first. Hence, I'm voting to close the question at this site. $\endgroup$
    – JimB
    Commented Jul 25, 2023 at 23:31
  • 1
    $\begingroup$ Also, I suspect (but could certainly be wrong) that there is little information added with the somewhat more vague knowledge that the PDF is monotonically decreasing. Your other question about adding knowledge of "radial symmetry" is different because that can be explicitly added to kernel density estimators. With a large sample size the estimated density will be close to being monotonically decreasing. $\endgroup$
    – JimB
    Commented Jul 25, 2023 at 23:39
  • $\begingroup$ You can do this in R with the Grenander Estimator: search.r-project.org/CRAN/refmans/fdrtool/html/grenander.html. $\endgroup$
    – JimB
    Commented Jul 25, 2023 at 23:50

4 Answers 4

7
$\begingroup$

If the Grenander Estimator is what you're looking for, then here is an implementation using Mathematica.

n = 100;  (* Sample size *)
SeedRandom[12345];
vals = RandomVariate[#, n] &@MarchenkoPasturDistribution[1];

(* Construct an empircal cdf *)
data = Transpose[{Sort[vals], ConstantArray[1, n]}];
data[[All, 2]] = Accumulate[data[[All, 2]]]/n;

(* Find the smallest concave majorant of the empirical cdf *)
(* This algorithm is probably fragile in that I'm assuming
a particular structure for the result of the ConvexHullRegion function.
Essentially I'm getting rid of the convex hull line segments "underneath"
the smallest concave majorant and including the point {0,0}. *)
ch = ConvexHullRegion[data];
max = Position[ch[[1, All, 2]], Max[ch[[1, All, 2]]]][[1, 1]];
min = Position[ch[[1, All, 2]], Min[ch[[1, All, 2]]]][[1, 1]];
ch = Join[ch[[1, max ;;]], {ch[[1, min]]}, {{0, 0}}];

(* Grenander CDF estimate *)
gCDF[x_] := Piecewise[Join[{{1, x >= Max[vals]}}, 
  Table[{ch[[i, 2]] + (ch[[i - 1, 2]] - ch[[i, 2]]) (x - ch[[i, 1]])/(ch[[i - 1, 1]] - ch[[i, 1]]), 
  ch[[i, 1]] <= x < ch[[i - 1, 1]]}, {i, 2, Length[ch]}]], 0]

(* Grenander PDF estimate *)
gPDF[x_] := Piecewise[Table[{(ch[[i, 2]] - ch[[i - 1, 2]])/(ch[[i, 1]] - ch[[i - 1, 1]]), 
  ch[[i, 1]] <= x <= ch[[i - 1, 1]]}, {i, 2, Length[ch]}], 0]

(* Plot pdfs *)
Plot[{gPDF[x], PDF[MarchenkoPasturDistribution[1], x],
  PDF[SmoothKernelDistribution[vals, Automatic, {"Bounded", {0, ∞}, "Gaussian"}], x]},
  {x, 0, 1.1 Max[vals]}, Exclusions -> None, PlotRange -> {All, {0, 2}},
  PlotStyle -> {Green, Red, Black}, Frame -> True,
  PlotLegends -> {"Grenander Estimate", "True density", "Smooth kernel density"},
  PlotLabel -> Style["PDF", 18, Bold, Black]]

PDF's

(* Plot cdfs *)
Plot[{gCDF[x], CDF[EmpiricalDistribution[vals], x], CDF[MarchenkoPasturDistribution[1], x]},
 {x, 0, 1.1 Max[vals]}, PlotLabel -> Style["CDF", 18, Bold, Black],
 PlotRange -> {{0, 1.1 Max[vals]}, {0, 1}}, Frame -> True, Exclusions -> None,
 PlotStyle -> {Green, Purple, Red}, PlotPoints -> 100,
 PlotLegends -> {"Grenander Estimate", "Empirical CDF", "True CDF"}]

CDF's

$\endgroup$
3
  • 2
    $\begingroup$ This is great! This was new to me so I did quite a bit of background study first; I found a lot of mentions of the Grenander Estimator in relevant publications and presentations, but could not find an implementation in mma. If you can work out the fragility kinks, then maybe this might be something worth putting on the wolfram function repository (it also seems this may be widely applicable to more problems)? $\endgroup$
    – ydd
    Commented Jul 27, 2023 at 1:10
  • $\begingroup$ @ydd Thanks. I'll have to think more about making the algorithm more stable. The aspects that do appear stable are (1) the convex hull is correct (not surprised at that) and (2) the points are in counter-clockwise order. $\endgroup$
    – JimB
    Commented Jul 27, 2023 at 18:28
  • $\begingroup$ @tdd Note: I now don't think my addition of the (0,0) point is appropriate. While that doesn't appear to matter too much, I'll check on the appropriateness of that shortly. I have figured out how to get the right "concave majorant of the empirical cdf" from Mathematica's ConvexHullRegion function while only relying on the assumption that the points are ordered in counter-clockwise order. I'll add that in the next few days. $\endgroup$
    – JimB
    Commented Jul 28, 2023 at 21:55
4
$\begingroup$

So I think I may be completely missing the mark here (if I am then I can delete this answer in a little bit) and this could be a much more complicated question than I'm qualified to answer, but you could start with the PDF of the HistogramDistribution:

SeedRandom[1];
vals = RandomVariate[#, 100] &@MarchenkoPasturDistribution[1];

hist = HistogramDistribution[vals, "Scott"]
DiscretePlot[PDF[hist, x], {x, 0, 4, 0.01}]

enter image description here

You can see the PDF is still not monotonically decreasing, but we can extract the PDF values at each bin, and the find the least squares solution that is monotonically decreasing:


pdfVals = hist["PDFValues"];
coefs = Array[c, Length@pdfVals];
constraints = Table[c[i + 1] <= c[i], {i, Length@coefs - 1}];
objectiveFunction = Norm[ pdfVals - coefs]^2;
forFit = Join[{objectiveFunction}, constraints];
fittedDecreasing = ArgMin[forFit, coefs];

delims = hist["BinDelimiters"];
booles = 
  Table[Boole[delims[[i]] <= x < delims[[i + 1]]], {i, 
    Length@delims - 1}];
fittedPDF = fittedDecreasing . booles;
DiscretePlot[{PDF[hist, x], fittedPDF}, {x, 0, 4, 0.01}]

enter image description here

$\endgroup$
7
  • $\begingroup$ Thats a neat trick, this might work with manual density estimation too (without binning, density is just the slope between successive points) $\endgroup$ Commented Jul 26, 2023 at 0:39
  • $\begingroup$ Ah so the CDF of the EmpiricalDistribution is just a piecewise function that increases by 1/100 at each value of vals right? So the slope at each point of the CDF is something like 0.01/Differences[Sort@vals] right? And you want the differences of the density to all be less than or equal to zero (monotonically decreasing), so I think this is the same as finding the best fit for 0.01/Differences[Sort@vals,2] (where Differences[_,2] is the second order difference) where everything is <= 0 right? Does that make sense? $\endgroup$
    – ydd
    Commented Jul 26, 2023 at 1:19
  • $\begingroup$ I modified the question to contain an example of computing empirical PDF without binning. I'm hoping monotonicity constraint can make the plot much better without introducing edge artifacts (like SmoothHistogram[vals]) $\endgroup$ Commented Jul 26, 2023 at 2:30
  • $\begingroup$ @YaroslavBulatov You can get rid of the edge artifacts with SmoothKernelDistribution[vals, Automatic, {"Bounded", {0, \[Infinity]}, "Gaussian"}] if what you mean by artifacts is the positive density to the left of zero. Until there is a large sample size, you won't have a monotone decreasing pdf. $\endgroup$
    – JimB
    Commented Jul 26, 2023 at 2:52
  • $\begingroup$ @JimB thanks for the pointer to the Grenander Estimator, this visualization seems to be what I need here $\endgroup$ Commented Jul 26, 2023 at 2:56
4
$\begingroup$
$Version

(* "13.3.0 for Mac OS X ARM (64-bit) (June 3, 2023)" *)

Clear["Global`*"]

SeedRandom[1234];

vals = RandomVariate[#, 100] &@MarchenkoPasturDistribution[1];

findDist = (# /. FindDistributionParameters[vals, #] &);

Pick several distributions that are monotonic for at least some values of their parameters.

distrSym = {
   MarchenkoPasturDistribution[λ, σ],
   MarchenkoPasturDistribution[λ],
   ChiSquareDistribution[ν],
   GammaDistribution[α, β],
   ExponentialDistribution[λ]};

Fit distributions to data

distr = findDist /@ distrSym

(* {MarchenkoPasturDistribution[0.980962, 1.01724], 
 MarchenkoPasturDistribution[0.986602, 1], ChiSquareDistribution[1.17914], 
 GammaDistribution[0.63868, 1.68378], ExponentialDistribution[0.929888]} *)

labels = StringReplace[ToString /@ distr, "Distribution" -> ""];

fit = NumberForm[#, {5, 3}] & /@ (DistributionFitTest[vals, #] & /@ distr)

enter image description here

colors = ColorData[97, "ColorList"];

Manipulate[
 d = Sort[d];
 Show[
  Histogram[vals, Automatic, "PDF"],
  Plot[
   Evaluate[Tooltip[PDF[#[[1]], x],
       StringForm["``, fit = ``", #[[1]], #[[2]]]] & /@
     Transpose[{distr, fit}][[d]]], {x, 0, 4},
   PlotStyle -> colors[[d]],
   PlotLegends -> Placed[
     LineLegend[
      StringForm["``, fit = ``", #[[1]], NumberForm[#[[2]], {5, 3}]] & /@ 
       Transpose[{labels, fit}][[d]]],
     {.6, .7}]],
  PlotRangeClipping -> False,
  ImageSize -> Large],
 {{d, Range[5], "distributions"},
  Thread[Range[Length[distr]] -> labels],
  ControlType -> TogglerBar,
  Appearance -> "Vertical" -> {Automatic, 3}}]

enter image description here

Select the distribution with the best fit.

$\endgroup$
3
$\begingroup$

I'd like to propose a radically simple (to the point of being stupid) sort of answer.

Why not use the uniform distribution as your kernel? Sounds crazy, but the UniformDistribution[{0, 2 v}] is monotonically decreasing and has a mean value of "v". When you sum them together you get a monotonically decreasing function. Trying it out:

SeedRandom[12345];
vals = RandomVariate[#, 100] &@MarchenkoPasturDistribution[1];

myEmpiricalPDF[x_] = 
 (vals /. v_Real -> PDF[UniformDistribution[{0, 2 v}], x] // Mean)

Okay, so yes that is an absurd definition and I wouldn't use it in this form in practice. We could write something much better, but if we just look at the result:

Plot[{myEmpiricalPDF[x], PDF[MarchenkoPasturDistribution[1], x]}, {x, 
  0, 5}]

enter image description here

Doesn't seem bad for the simplest possible kernel (although it wouldn't actually converge to the correct distribution in this case)

I had tons of fun thinking about this, so thanks for posting it. I considered using several other kernels that would give smooth results with the same conceptual plan, but this seems to work okay.

$\endgroup$
3
  • 1
    $\begingroup$ Wow, that plot looks pretty nice. My original motivation was that empirical CDF looks nice, but empirical PDF does not. I have to think some more why this trick works.... $\endgroup$ Commented Jul 28, 2023 at 17:58
  • 1
    $\begingroup$ Interesting. This spreads out the influence of larger observations proportionally more than smaller observations somewhat analogous to an adaptive kernel estimator. But there doesn't appear to be any control for how much spread is allowed. I've tried your procedure with larger sample sizes but the fit appears worse as the sample size increases. But no matter how "crazy" (your words not mine) it appears, the long-term characteristics are what matter. Have you looked into long-term characteristics of this estimator? $\endgroup$
    – JimB
    Commented Jul 28, 2023 at 21:52
  • $\begingroup$ @jimb I had the same thought, but after sleeping on it, I've decided the bug is actually a feature. We lose resolution on larger values proportional to their value and it's sorta appropriate to lose resolution in the tail of the distribution. $\endgroup$
    – Searke
    Commented Jul 29, 2023 at 20:54

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.