text: #!/usr/bin/env python3
  
  """
  Code for "Python Histogram Plotting: NumPy, Matplotlib, Pandas & Seaborn".
  
  https://realpython.com/python-histograms/
  
  To run this script interactively:
  $ ipython ./histograms.py -i
  or
  $ python3 -i ./histograms.py
  """
  
  from collections import Counter
  import random
  import sys
  import warnings
  
  import matplotlib.pyplot as plt
  import numpy as np
  import pandas as pd
  from scipy import stats
  import seaborn as sns
  
  print(__doc__)
  
  
  random.seed(1)
  np.random.seed(444)
  np.set_printoptions(precision=3)
  plt.ioff()
  
  # ---------------------------------------------------------------------
  # Histograms in Pure Python
  # ---------------------------------------------------------------------
  
  
  def count_elements(seq) -> dict:
      """Tally elements from `seq`."""
      hist = {}
      # Bind method to a variable for faster calls within loop.
      get = hist.get
      for i in seq:
          hist[i] = get(i, 0) + 1
      return hist
  
  
  a = (0, 1, 1, 1, 2, 3, 7, 7, 23)
  counted = count_elements(a)
  recounted = Counter(a)
  # Both `.items()` are instances of Python builtin `dict_items`.
  assert counted.items() == recounted.items(), "Frequency tables unequal."
  print("`a`:", a)
  print("Manual frequency table of `a`:", counted)
  print("With `collections.Counter`:", recounted)
  print()
  
  
  def ascii_histogram(seq) -> None:
      """A horizontal frequency-table/histogram plot."""
      counted = count_elements(seq)
      for k in sorted(counted):
          print("{0:5d} {1}".format(k, "+" * counted[k]))
  
  
  vals = [1, 3, 4, 6, 8, 9, 10]
  # Each number in `vals` will occur between 5 and 15 times.
  freq = [random.randint(5, 15) for _ in vals]
  
  data = []
  for f, v in zip(freq, vals):
      data.extend([v] * f)
  print("ASCII histogram of `data`:")
  ascii_histogram(data)
  print()
  
  
  # ---------------------------------------------------------------------
  # Histograms Calculations in NumPy
  # ---------------------------------------------------------------------
  
  d = np.random.laplace(loc=15, scale=3, size=500)
  hist, bin_edges = np.histogram(d)
  print("first 5 elements of `d`:", d[:5])
  print("hist:", hist)
  print("bin_edges:", bin_edges)
  
  bcounts = np.bincount(a)
  hist, _ = np.histogram(a, range=(0, max(a)), bins=max(a) + 1)
  print(bcounts)
  assert np.array_equal(hist, bcounts), "Bincounts unequal."
  
  # Reproducing `collections.Counter`
  print(
      "Reproducing `collections.Counter`:",
      dict(zip(np.unique(a), bcounts[bcounts.nonzero()])),
  )
  
  
  # ---------------------------------------------------------------------
  # Visualizing Histograms with Matplotlib & Pandas
  # ---------------------------------------------------------------------
  
  n, bins, patches = plt.hist(
      x=d, bins="auto", color="#0504aa", alpha=0.7, rwidth=0.85
  )
  plt.grid(axis="y", alpha=0.75)
  plt.xlabel("Value")
  plt.ylabel("Frequency")
  plt.title("My Very Own Histogram")
  plt.text(23, 45, r"$\mu=15, b=3$")
  maxfreq = n.max()
  # Set a clean upper y-axis limit.
  plt.ylim(ymax=np.ceil(maxfreq / 10) * 10 if maxfreq % 10 else maxfreq + 10)
  plt.show()
  
  
  # Generate data on commute times.
  size, scale = 1000, 10
  commutes = pd.Series(np.random.gamma(scale, size=size) ** 1.5)
  
  commutes.plot.hist(bins=15, rwidth=0.9, color="#607c8e")
  plt.title("Commute Times for 1,000 Commuters")
  plt.xlabel("Commute Time (min)")
  plt.ylabel("Counts")
  plt.grid(axis="y", alpha=0.75)
  plt.show()
  
  # ---------------------------------------------------------------------
  # Plotting a Kernel Density Estimate
  # ---------------------------------------------------------------------
  
  # Sample from two different normal distributions.
  means = 10, 20
  stdevs = 4, 2
  dist = pd.DataFrame(
      np.random.normal(loc=means, scale=stdevs, size=(1000, 2)),
      columns=["a", "b"],
  )
  print(dist.agg(["min", "max", "mean", "std"]).round(decimals=2))
  
  fig, ax = plt.subplots()
  dist.plot.kde(ax=ax, grid=True, legend=False, title="Histogram: A vs. B")
  dist.plot.hist(density=True, ax=ax)
  ax.set_ylabel("Probability")
  ax.set_facecolor("#d8dcd6")
  ax.grid(axis="y")
  plt.show()
  
  
  # An object representing the "frozen" analytical distribution.
  # Defaults to the standard normal distribution, N~(0, 1)
  dist = stats.norm()
  
  # Draw random samples from the population you built above.
  # This is just a sample, so the mean and std. deviation should
  # be close to (1, 0).
  samp = dist.rvs(size=1000)
  
  # `ppf()`: percent point function (inverse of cdf — percentiles).
  x = np.linspace(start=stats.norm.ppf(0.01), stop=stats.norm.ppf(0.99), num=250)
  gkde = stats.gaussian_kde(dataset=samp)
  
  # `gkde.evaluate()` estimates the PDF itself.
  fig, ax = plt.subplots()
  ax.plot(
      x,
      dist.pdf(x),
      linestyle="solid",
      c="red",
      lw=3,
      alpha=0.8,
      label="Analytical (True) PDF",
  )
  ax.plot(
      x,
      gkde.evaluate(x),
      linestyle="dashed",
      c="black",
      lw=2,
      label="PDF Estimated via KDE",
  )
  ax.legend(loc="best", frameon=False)
  ax.set_title("Analytical vs. Estimated PDF")
  ax.set_ylabel("Probability")
  ax.text(-2.0, 0.35, r"$f(x) = \frac{\exp(-x^2/2)}{\sqrt{2*\pi}}$", fontsize=12)
  plt.show()
  
  # ---------------------------------------------------------------------
  # A Fancy Alternative with Seaborn
  # ---------------------------------------------------------------------
  
  sns.set_style("darkgrid")
  # Suppress the kwarg warning related to normed/density from Matplotlib.
  with warnings.catch_warnings():
      warnings.simplefilter("ignore", category=UserWarning)
  
      sns.distplot(d)
      plt.title("Seaborn's distplot()")
      plt.show()
  
      sns.distplot(d, fit=stats.laplace, kde=False)
      plt.title("Histogram with Fitted Laplace Distribution")
      plt.show()
  
  
  data = np.random.choice(
      np.arange(10), size=10000, p=np.linspace(1, 11, 10) / 60
  )
  s = pd.Series(data)
  print(s.value_counts())
  print(s.value_counts(normalize=True).head())
  
  ages = pd.Series(
      [1, 1, 3, 5, 8, 10, 12, 15, 18, 18, 19, 20, 25, 30, 40, 51, 52]
  )
  bins = (0, 10, 13, 18, 21, np.inf)  # The edges
  labels = ("child", "preteen", "teen", "military_age", "adult")
  groups, _ = pd.cut(ages, bins=bins, labels=labels, retbins=True)
  print(groups.value_counts())
  print(pd.concat((ages, groups), axis=1).rename(columns={0: "age", 1: "group"}))
  
  sys.exit()
  