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

ClearAll[hangingRootogram]
hangingRootogram[dat_, estdist_, binspec_: Automatic][sc___ : .9, o: OptionsPattern[]] :=
With[{hd = HistogramDistribution[dat, binspec], bins = HistogramList[dat, binspec][]},
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, 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 '16 at 15:01
• So much from a compact function. One of the coolest answers I've seen in a while. +1
– ciao
Jul 8 '16 at 7:22
• Thank you @ciao for the kind words and the upvote.
– kglr
Jul 8 '16 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[], Line[p_] :> p, \[Infinity]][], 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[], 2, 1], ps[]}];
If[sqRoot, ps[[All, 2]] = Sqrt[ps[[All, 2]]]];
s = Max[Abs[Differences[ps[[All, 1]]]]];
Graphics[{
GrayLevel[0.7],
Map[Rectangle[{#[] - s/2.5,
shpoints[[nf[#[]][], 2]] - #[]}, {#[] + s/2.5,
shpoints[[nf[#[]][], 2]]}] &, ps],
Blue, Line[Select[shpoints, x0 <= #[] <= x1 &]],
Red, Point[Map[shpoints[[nf[#[]][]]] &, ps]]}, opts,
Axes -> True, AspectRatio -> 1/GoldenRatio]
];

dist = PoissonDistribution;
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 '16 at 13:19
• @J.M. Yeah, it would be nice to read that article... Jul 7 '16 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 '16 at 15:19
• @Szabolcs related to your comment: (8652), (16694), (49331), (111438) Jul 8 '16 at 6:39
• Aha, found it! It wasn't in EDA after all, upon checking my copy. Jul 10 '16 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];
{
Opacity@.5,
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 '16 at 14:13
• @AntonAntonov Thanks for attention, I will try to investigate later, have to go now.
– Kuba
Jul 7 '16 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 '16 at 15:04
• @Szabolcs Please see my answer -- it has graphs/histograms with PoissonDistribution. Jul 7 '16 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 '16 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] 