# How to draw a hanging rootogram in Mathematica?

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.

• 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. Jul 7, 2016 at 11:58
• "some data" - can you please show us this? Jul 7, 2016 at 12:40

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}]]]}]


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]}]


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}]]]}]


• 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. Jul 7, 2016 at 15:01
• So much from a compact function. One of the coolest answers I've seen in a while. +1
– ciao
Jul 8, 2016 at 7:22
• Thank you @ciao for the kind words and the upvote.
– kglr
Jul 8, 2016 at 8:44

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]}}]


• I think the best way to answer this would be to look up Tukey's description in EDA. I'd check myself, but I'm far away from my copy. :( Jul 7, 2016 at 13:19
• @J.M. Yeah, it would be nice to read that article... Jul 7, 2016 at 14:48
• I noticed {_?NumberQ..} here. You may be interested in this chat from yesterday: chat.stackexchange.com/transcript/message/30840228#30840228 Jul 7, 2016 at 15:19
• @Szabolcs related to your comment: (8652), (16694), (49331), (111438) Jul 8, 2016 at 6:39
• Aha, found it! It wasn't in EDA after all, upon checking my copy. Jul 10, 2016 at 21:31

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
]


Of course the more points the better match:

n = 10000


• I think this does not work with well with other distributions like PoissonDistribution or WeibullDistribution. Jul 7, 2016 at 14:13
• @AntonAntonov Thanks for attention, I will try to investigate later, have to go now.
– Kuba
Jul 7, 2016 at 14:25
• @AntonAntonov Seems to work well with Weibull. Poisson is a discrete distribution so neither Plot nor a non-integer binned histogram would be appropriate. Jul 7, 2016 at 15:04
• @Szabolcs Please see my answer -- it has graphs/histograms with PoissonDistribution. Jul 7, 2016 at 15:10
• @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. Jul 7, 2016 at 15:18

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]