text: r"""
  Jackknife Calculations of Error on Mean
  ---------------------------------------
  Figure 4.4.
  
  The jackknife uncertainty estimates for the width of a Gaussian distribution.
  This example uses the same data as figure 4.3. The upper panel shows a
  histogram of the widths determined using the sample standard deviation, and
  using the interquartile range. The lower panel shows the corrected jackknife
  estimates (eqs. 4.33 and 4.35) for the two methods. The gray lines show the
  theoretical results, given by eq. 3.35 for :math:`\sigma` and eq. 3.37 for
  :math:`\sigma_G`. The result for :math:`\sigma` matches the theoretical
  result almost exactly, but note the failure of the jackknife to correctly
  estimate :math:`\sigma_G` (see the text for a discussion of this result).
  """
  # Author: Jake VanderPlas
  # License: BSD
  #   The figure produced by this code is published in the textbook
  #   "Statistics, Data Mining, and Machine Learning in Astronomy" (2013)
  #   For more information, see http://astroML.github.com
  #   To report a bug or issue, use the following forum:
  #    https://groups.google.com/forum/#!forum/astroml-general
  import numpy as np
  from scipy.stats import norm
  from matplotlib import pyplot as plt
  
  #----------------------------------------------------------------------
  # This function adjusts matplotlib settings for a uniform feel in the textbook.
  # Note that with usetex=True, fonts are rendered with LaTeX.  This may
  # result in an error if LaTeX is not installed on your system.  In that case,
  # you can set usetex to False.
  from astroML.plotting import setup_text_plots
  setup_text_plots(fontsize=8, usetex=True)
  
  #------------------------------------------------------------
  # sample values from a normal distribution
  np.random.seed(123)
  m = 1000  # number of points
  data = norm(0, 1).rvs(m)
  
  #------------------------------------------------------------
  # Compute jackknife resamplings of data
  from astroML.resample import jackknife
  from astroML.stats import sigmaG
  
  # mu1 is the mean of the standard-deviation-based width
  mu1, sigma_mu1, mu1_raw = jackknife(data, np.std,
                                      kwargs=dict(axis=1, ddof=1),
                                      return_raw_distribution=True)
  
  pdf1_theory = norm(1, 1. / np.sqrt(2 * (m - 1)))
  pdf1_jackknife = norm(mu1, sigma_mu1)
  
  # mu2 is the mean of the interquartile-based width
  #  WARNING: do not use the following in practice.  This example
  #           shows that jackknife fails for rank-based statistics.
  mu2, sigma_mu2, mu2_raw = jackknife(data, sigmaG,
                                      kwargs=dict(axis=1),
                                      return_raw_distribution=True)
  pdf2_theory = norm(data.std(), 1.06 / np.sqrt(m))
  pdf2_jackknife = norm(mu2, sigma_mu2)
  print mu2, sigma_mu2
  
  #------------------------------------------------------------
  # plot the results
  print "mu_1 mean: %.2f +- %.2f" % (mu1, sigma_mu1)
  print "mu_2 mean: %.2f +- %.2f" % (mu2, sigma_mu2)
  
  fig = plt.figure(figsize=(5, 2))
  fig.subplots_adjust(left=0.11, right=0.95, bottom=0.2, top=0.9,
                      wspace=0.25)
  
  ax = fig.add_subplot(121)
  ax.hist(mu1_raw, np.linspace(0.996, 1.008, 100),
          label=r'$\sigma^*\ {\rm (std.\ dev.)}$',
          histtype='stepfilled', fc='white', normed=False)
  ax.hist(mu2_raw, np.linspace(0.996, 1.008, 100),
          label=r'$\sigma_G^*\ {\rm (quartile)}$',
          histtype='stepfilled', fc='gray', normed=False)
  ax.legend(loc='upper left', handlelength=2)
  
  ax.xaxis.set_major_locator(plt.MultipleLocator(0.004))
  ax.set_xlabel(r'$\sigma^*$')
  ax.set_ylabel(r'$N(\sigma^*)$')
  ax.set_xlim(0.998, 1.008)
  ax.set_ylim(0, 550)
  
  ax = fig.add_subplot(122)
  x = np.linspace(0.45, 1.15, 1000)
  ax.plot(x, pdf1_jackknife.pdf(x),
          color='blue', ls='dashed', label=r'$\sigma\ {\rm (std.\ dev.)}$',
          zorder=2)
  ax.plot(x, pdf1_theory.pdf(x), color='gray', zorder=1)
  ax.plot(x, pdf2_jackknife.pdf(x),
          color='red', label=r'$\sigma_G\ {\rm (quartile)}$', zorder=2)
  ax.plot(x, pdf2_theory.pdf(x), color='gray', zorder=1)
  plt.legend(loc='upper left', handlelength=2)
  
  ax.set_xlabel(r'$\sigma$')
  ax.set_ylabel(r'$p(\sigma|x,I)$')
  ax.set_xlim(0.45, 1.15)
  ax.set_ylim(0, 24)
  
  plt.show()
  