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
    $\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$ – Anton Antonov Jul 7 '16 at 11:58
  • 2
    $\begingroup$ "some data" - can you please show us this? $\endgroup$ – J. M. will be back soon Jul 7 '16 at 12:40
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$
  • $\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$ – Anton Antonov Jul 7 '16 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 '16 at 7:22
  • $\begingroup$ Thank you @ciao for the kind words and the upvote. $\endgroup$ – kglr Jul 8 '16 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$
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];
     {
       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
]

enter image description here

Of course the more points the better match:

n = 10000

enter image description here

$\endgroup$
  • $\begingroup$ I think this does not work with well with other distributions like PoissonDistribution or WeibullDistribution. $\endgroup$ – Anton Antonov Jul 7 '16 at 14:13
  • $\begingroup$ @AntonAntonov Thanks for attention, I will try to investigate later, have to go now. $\endgroup$ – Kuba Jul 7 '16 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 '16 at 15:04
  • $\begingroup$ @Szabolcs Please see my answer -- it has graphs/histograms with PoissonDistribution. $\endgroup$ – Anton Antonov Jul 7 '16 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 '16 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, privacy policy and cookie policy

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