code: """Stack loss data"""
  
  __all__ = ['COPYRIGHT','TITLE','SOURCE','DESCRSHORT','DESCRLONG','NOTE', 'load']
  
  
  __docformat__ = 'restructuredtext'
  
  COPYRIGHT   = """This is public domain. """
  TITLE       = __doc__
  SOURCE      = """
  Brownlee, K. A. (1965), "Statistical Theory and Methodology in
  Science and Engineering", 2nd edition, New York:Wiley.
  """
  
  DESCRSHORT  = """Stack loss plant data of Brownlee (1965)"""
  
  DESCRLONG   = """The stack loss plant data of Brownlee (1965) contains
  21 days of measurements from a plant's oxidation of ammonia to nitric acid.
  The nitric oxide pollutants are captured in an absorption tower."""
  
  NOTE        = """
  Number of Observations - 21
  
  Number of Variables - 4
  
  Variable name definitions::
  
      STACKLOSS - 10 times the percentage of ammonia going into the plant that
                  escapes from the absoroption column
      AIRFLOW   - Rate of operation of the plant
      WATERTEMP - Cooling water temperature in the absorption tower
      ACIDCONC  - Acid concentration of circulating acid minus 50 times 10.
  """
  
  from numpy import recfromtxt, column_stack, array
  from scikits.statsmodels.datasets import Dataset
  from os.path import dirname, abspath
  
  def load():
      """
      Load the stack loss data and returns a Dataset class instance.
  
      Returns
      --------
      Dataset instance:
          See DATASET_PROPOSAL.txt for more information.
      """
      filepath = dirname(abspath(__file__))
      data = recfromtxt(open(filepath + '/stackloss.csv',"rb"), delimiter=",",
              names=True, dtype=float)
      names = list(data.dtype.names)
      endog = array(data[names[0]], dtype=float)
      endog_name = names[0]
      exog = column_stack(data[i] for i in names[1:]).astype(float)
      exog_name = names[1:]
      dataset = Dataset(data=data, names=names, endog=endog, exog=exog,
              endog_name = endog_name, exog_name=exog_name)
      return dataset
   import numpy as np
  from scipy.stats import scoreatpercentile as sap
  
  #from scipy.stats import norm
  
  def _select_sigma(X):
      """
      Returns the smaller of std(X, ddof=1) or normalized IQR(X) over axis 0.
  
      References
      ----------
      Silverman (1986) p.47
      """
  #    normalize = norm.ppf(.75) - norm.ppf(.25)
      normalize = 1.349
  #    IQR = np.subtract.reduce(percentile(X, [75,25],
  #                             axis=axis), axis=axis)/normalize
      IQR = (sap(X, 75) - sap(X, 25))/normalize
      return np.minimum(np.std(X, axis=0, ddof=1), IQR)
  
  
  ## Univariate Rule of Thumb Bandwidths ##
  def bw_scott(x):
      """
      Scott's Rule of Thumb
  
      Parameter
      ---------
      x : array-like
          Array for which to get the bandwidth
  
      Returns
      -------
      bw : float
          The estimate of the bandwidth
  
      Notes
      -----
      Returns 1.059 * A * n ** (-1/5.)
  
      A = min(std(x, ddof=1), IQR/1.349)
      IQR = np.subtract.reduce(np.percentile(x, [75,25]))
  
      References
      ---------- ::
  
      Scott, D.W. (1992) `Multivariate Density Estimation: Theory, Practice, and
          Visualization.`
      """
      A = _select_sigma(x)
      n = len(x)
      return 1.059 * A * n ** -.2
  
  def bw_silverman(x):
      """f
      Silverman's Rule of Thumb
  
      Parameter
      ---------
      x : array-like
          Array for which to get the bandwidth
  
      Returns
      -------
      bw : float
          The estimate of the bandwidth
  
      Notes
      -----
      Returns .9 * A * n ** (-1/5.)
  
      A = min(std(x, ddof=1), IQR/1.349)
      IQR = np.subtract.reduce(np.percentile(x, [75,25]))
  
      References
      ---------- ::
  
      Silverman, B.W. (1986) `Density Estimation.`
      """
      A = _select_sigma(x)
      n = len(x)
      return .9 * A * n ** -.2
  
  ## Plug-In Methods ##
  
  ## Least Squares Cross-Validation ##
  
  ## Helper Functions ##
  
  bandwidth_funcs = dict(scott=bw_scott,silverman=bw_silverman)
  
  def select_bandwidth(X, bw, kernel):
      """
      Selects bandwidth
      """
      bw = bw.lower()
      if bw not in ["scott","silverman"]:
          raise ValueError("Bandwidth %s not understood" % bw)
  #TODO: uncomment checks when we have non-rule of thumb bandwidths for diff. kernels
  #    if kernel == "gauss":
      return bandwidth_funcs[bw](X)
  #    else:
  #        raise ValueError("Only Gaussian Kernels are currently supported")
  
  
apis: numpy.array numpy.column_stack
apis: scipy.stats.scoreatpercentile numpy.std