8
$\begingroup$

I am trying to plot a hanging rootogram of some data in Mathematica. I can't seem to find a built in function for it, while simply using Histogram (on "transformed" data) does not seem to plot what I want.

$\endgroup$
2
  • 2
    $\begingroup$ By "a hanging rootogram" do you mean like the one shown here ? If not, please change the link I put in your question to something more relevant. $\endgroup$ Jul 7, 2016 at 11:58
  • 2
    $\begingroup$ "some data" - can you please show us this? $\endgroup$ Jul 7, 2016 at 12:40

4 Answers 4

18
$\begingroup$
ClearAll[hangingRootogram]
hangingRootogram[dat_, estdist_, binspec_: Automatic][sc___ : .9, o: OptionsPattern[]] := 
 With[{hd = HistogramDistribution[dat, binspec], bins = HistogramList[dat, binspec][[1]]}, 
  With[{es = sc  Min@Differences@bins}, 
   DiscretePlot[{Sqrt@PDF[estdist, x] - Sqrt@PDF[hd, x], Sqrt@PDF[estdist, x]}, {x, bins}, 
    ExtentSize -> es, PlotMarkers -> {None, {"Point", Large}}, Joined -> {False, True}, 
    Filling -> {1 -> {2}}, o]]]

Examples:

data = RandomVariate[NegativeBinomialDistribution[10, 0.3], 10^2];
edist = EstimatedDistribution[data, NegativeBinomialDistribution[n, p], 
  ParameterEstimator -> "MethodOfMoments"];

Row[{Histogram[data, Automatic, "PDF", ImageSize -> 400], 
  hangingRootogram[data, edist][.8, ImageSize -> 400, 
   PlotStyle -> {Blue, Red}, 
   FillingStyle -> Directive[Opacity[.7], Blue, EdgeForm[{Blue, Thick}]]]}]

Mathematica graphics

Row[hangingRootogram[data, EstimatedDistribution[data, #, 
      ParameterEstimator -> "MethodOfMoments"]][.8, ImageSize -> 400, 
    PlotStyle -> {Blue, Red}, PlotLabel -> #, 
    FillingStyle -> Directive[Opacity[.7], Blue, EdgeForm[{Blue, Thick}]], 
    PlotRange -> Full] & /@ 
   {NegativeBinomialDistribution[n, p], NegativeBinomialDistribution[n, .5],
    PoissonDistribution[n]}]

Mathematica graphics

data = RandomVariate[PoissonDistribution[5], 10^3];
edist = EstimatedDistribution[data, PoissonDistribution[n], 
  ParameterEstimator -> "MethodOfMoments"];

Row[{Histogram[data, {3}, "PDF", ImageSize -> 400], 
  hangingRootogram[data, edist, {3}][.8, ImageSize -> 400, 
   PlotStyle -> {Blue, Red}, 
   FillingStyle -> Directive[Opacity[.7], Blue, EdgeForm[{Blue, Thick}]]]}]

Mathematica graphics

$\endgroup$
3
  • $\begingroup$ Probably the best solution so far (+1). I think it complements mine since it uses PDFs and use bin counts. Please consider overloading hangingRootogram to not take the distribution parameter edist. $\endgroup$ Jul 7, 2016 at 15:01
  • 2
    $\begingroup$ So much from a compact function. One of the coolest answers I've seen in a while. +1 $\endgroup$
    – ciao
    Jul 8, 2016 at 7:22
  • $\begingroup$ Thank you @ciao for the kind words and the upvote. $\endgroup$
    – kglr
    Jul 8, 2016 at 8:44
9
$\begingroup$

The defined function RootHistogram makes a "hanging rootogram" more-or-less following this definition.

The first argument is the data. The second argument dist is optional distribution. The function uses SmoothHistogram for the hanging curve and the third argument, bandWidth, is the band width argument of SmoothHistogram. The bspec argument is given to HistogramList. The sqRoot argument is in adherence to the mentioned definition:

[...] As in the rootogram, the vertical axis is scaled to the square-root of the frequencies so as to draw attention to discrepancies in the tails of the distribution.

Clear[RootHistogram]
RootHistogram[data : {_?NumberQ ..}, dist_: Automatic, 
   bandWidth_: "StandardDeviation", bspec_: Automatic, 
   sqRoot : (True | False) : True, opts : OptionsPattern[]] :=
  Block[{gr, shpoints, nf, x0, x1, s, xs, ds, ps},
   gr = SmoothHistogram[data, bandWidth, "Intensity"];
   shpoints = 
    SortBy[Cases[gr[[1]], Line[p_] :> p, \[Infinity]][[1]], First];
   If[! TrueQ[dist === Automatic],
    ds = Table[PDF[dist, x], {x, shpoints[[All, 1]]}];
    ds = Rescale[ds, MinMax[ds], MinMax[shpoints[[All, 2]]]];
    shpoints[[All, 2]] = ds
    ];
   If[sqRoot, shpoints[[All, 2]] = Sqrt[shpoints[[All, 2]]]];
   nf = Nearest[shpoints[[All, 1]] -> Automatic];
   {x0, x1} = MinMax[data];
   ps = HistogramList[data, bspec];
   ps = Transpose[{Mean /@ Partition[ps[[1]], 2, 1], ps[[2]]}];
   If[sqRoot, ps[[All, 2]] = Sqrt[ps[[All, 2]]]];
   s = Max[Abs[Differences[ps[[All, 1]]]]];
   Graphics[{
     GrayLevel[0.7],
     Map[Rectangle[{#[[1]] - s/2.5, 
         shpoints[[nf[#[[1]]][[1]], 2]] - #[[2]]}, {#[[1]] + s/2.5, 
         shpoints[[nf[#[[1]]][[1]], 2]]}] &, ps],
     Blue, Line[Select[shpoints, x0 <= #[[1]] <= x1 &]],
     Red, Point[Map[shpoints[[nf[#[[1]]][[1]]]] &, ps]]}, opts, 
    Axes -> True, AspectRatio -> 1/GoldenRatio]
   ];

dist = PoissonDistribution[8];
data = RandomVariate[dist, 500];

opts = {ImageSize -> 450, Axes -> False, Frame -> True};
Grid[{{Histogram[data, 20, PlotLabel -> "Histogram", opts],
   RootHistogram[data, Automatic, "StandardDeviation", 20, True, 
    PlotLabel -> 
     "\!\(\*SqrtBox[\(SmoothHistogram\)]\) with hanging \
\!\(\*SqrtBox[\(HistogramList\)]\) panels", opts]},
  {RootHistogram[data, Automatic, "StandardDeviation", 20, False, 
    PlotLabel -> "SmoothHistogram with hanging HistogramList panels", opts],
   RootHistogram[data, NormalDistribution[11, 2],"StandardDeviation", 20, True, PlotLabel -> 
     "\!\(\*SqrtBox[\(Max[SmoothHistogram] PDF[N[11, 2], x]\)]\) with \
hanging \!\(\*SqrtBox[\(HistogramList\)]\) panels", opts]}}]

enter image description here

$\endgroup$
8
7
$\begingroup$

I don't know how to interprete scaling of frequencies and associated expected curve so I will just plot PDF. This answer isn't complete then!

Here is a simple way to hang those bars using ChartElementFunction:

d = NormalDistribution[0, 1]
n = 100
data = RandomVariate[d, n];
bspec = {-5, 5, .5};

 f[{{xmin_, xmax_}, {ymin_, ymax_}}, ___] :=  Module[{
     m = Mean@{xmin, xmax}, yMax
     }, 
     yMax = PDF[d, m];
     {
       [email protected], 
       Translate[Rectangle[{xmin, ymin}, {xmax, ymax}], {0, yMax - ymax}],
       AbsolutePointSize@7, Red, 
       Point[{m, yMax}]
     }
];


Show[
 Plot[ PDF[d, x], {x, #, #2}] & @@ bspec,(*expected*)
 Histogram[data, bspec, (*experimental*) 
    "PDF", 
    ChartElementFunction -> f
 ],
 PlotRange -> All, 
 Frame -> True, 
 GridLines -> {{}, {0}}, 
 GridLinesStyle -> Thick
]

enter image description here

Of course the more points the better match:

n = 10000

enter image description here

$\endgroup$
7
  • $\begingroup$ I think this does not work with well with other distributions like PoissonDistribution or WeibullDistribution. $\endgroup$ Jul 7, 2016 at 14:13
  • $\begingroup$ @AntonAntonov Thanks for attention, I will try to investigate later, have to go now. $\endgroup$
    – Kuba
    Jul 7, 2016 at 14:25
  • $\begingroup$ @AntonAntonov Seems to work well with Weibull. Poisson is a discrete distribution so neither Plot nor a non-integer binned histogram would be appropriate. $\endgroup$
    – Szabolcs
    Jul 7, 2016 at 15:04
  • $\begingroup$ @Szabolcs Please see my answer -- it has graphs/histograms with PoissonDistribution. $\endgroup$ Jul 7, 2016 at 15:10
  • $\begingroup$ @Anton Yes, but yours simply does something different. It uses SmoothHistogram. Kuba uses PDF. That would require special handling for discrete distributions (which don't technically have a probability density). I don't think the two approaches are comparable in a work / doesn't work way. They do different things. $\endgroup$
    – Szabolcs
    Jul 7, 2016 at 15:18
2
$\begingroup$

This is a much simpler approach than already given and simply takes theoretical and measured values:

rootogram[theory_, observations_] := Show[{
   ListLinePlot[{theory}, PlotMarkers -> {Automatic, 10}],
   Graphics[{Table[
      Line[{{i, theory[[i]]}, {i, 
         measurements[[i]] - theory[[i]]}}], {i, Length[theory]}]}]
   }]

theory = {3, 5, 7, 9, 11, 13, 15, 17, 19, 21};
measurements = {2, 4, 7, 10, 12, 10, 16, 18, 19, 20};
rootogram[theory, measurements]

enter image description here

$\endgroup$

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.