title: AttributeError: module 'statsmodels' has no attribute 'datasets'
body: <p>I am using the code in this post:</p>
  
  <p><a href="https://stackoverflow.com/questions/6620471/fitting-empirical-distribution-to-theoretical-ones-with-scipy-python">Fitting empirical distribution to theoretical ones with Scipy (Python)?</a></p>
  
  <p>However, it keeps giving me this error:</p>
  
  <blockquote>
    <p>File "fitdistr.py", line 102, in 
        data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())
    AttributeError: module 'statsmodels' has no attribute 'datasets'</p>
  </blockquote>
  
  <p>Is there a way to load from a CSV file, as there was a comment to use this:</p>
  
  <p><a href="https://drive.google.com/file/d/1kOa35v0NHjWHK8VmMDY-sKiGZyeJzw6z/view" rel="nofollow noreferrer">https://drive.google.com/file/d/1kOa35v0NHjWHK8VmMDY-sKiGZyeJzw6z/view</a>
  Which is CSV data.</p>
  
  <pre><code>import warnings
  import numpy as np
  import pandas as pd
  import scipy.stats as st
  import statsmodels as sm
  import matplotlib
  import matplotlib.pyplot as plt
  
  matplotlib.rcParams['figure.figsize'] = (16.0, 12.0)
  matplotlib.style.use('ggplot')
  
  # Create models from data
  def best_fit_distribution(data, bins=200, ax=None):
      """Model data by finding best fit distribution to data"""
      # Get histogram of original data
      y, x = np.histogram(data, bins=bins, density=True)
      x = (x + np.roll(x, -1))[:-1] / 2.0
  
      # Distributions to check
      DISTRIBUTIONS = [        
          st.alpha,st.anglit,st.arcsine,st.beta,st.betaprime,st.bradford,st.burr,st.cauchy,st.chi,st.chi2,st.cosine,
          st.dgamma,st.dweibull,st.erlang,st.expon,st.exponnorm,st.exponweib,st.exponpow,st.f,st.fatiguelife,st.fisk,
          st.foldcauchy,st.foldnorm,st.frechet_r,st.frechet_l,st.genlogistic,st.genpareto,st.gennorm,st.genexpon,
          st.genextreme,st.gausshyper,st.gamma,st.gengamma,st.genhalflogistic,st.gilbrat,st.gompertz,st.gumbel_r,
          st.gumbel_l,st.halfcauchy,st.halflogistic,st.halfnorm,st.halfgennorm,st.hypsecant,st.invgamma,st.invgauss,
          st.invweibull,st.johnsonsb,st.johnsonsu,st.ksone,st.kstwobign,st.laplace,st.levy,st.levy_l,st.levy_stable,
          st.logistic,st.loggamma,st.loglaplace,st.lognorm,st.lomax,st.maxwell,st.mielke,st.nakagami,st.ncx2,st.ncf,
          st.nct,st.norm,st.pareto,st.pearson3,st.powerlaw,st.powerlognorm,st.powernorm,st.rdist,st.reciprocal,
          st.rayleigh,st.rice,st.recipinvgauss,st.semicircular,st.t,st.triang,st.truncexpon,st.truncnorm,st.tukeylambda,
          st.uniform,st.vonmises,st.vonmises_line,st.wald,st.weibull_min,st.weibull_max,st.wrapcauchy
      ]
  
      # Best holders
      best_distribution = st.norm
      best_params = (0.0, 1.0)
      best_sse = np.inf
  
      # Estimate distribution parameters from data
      for distribution in DISTRIBUTIONS:
  
          # Try to fit the distribution
          try:
              # Ignore warnings from data that can't be fit
              with warnings.catch_warnings():
                  warnings.filterwarnings('ignore')
  
                  # fit dist to data
                  params = distribution.fit(data)
  
                  # Separate parts of parameters
                  arg = params[:-2]
                  loc = params[-2]
                  scale = params[-1]
  
                  # Calculate fitted PDF and error with fit in distribution
                  pdf = distribution.pdf(x, loc=loc, scale=scale, *arg)
                  sse = np.sum(np.power(y - pdf, 2.0))
  
                  # if axis pass in add to plot
                  try:
                      if ax:
                          pd.Series(pdf, x).plot(ax=ax)
                      end
                  except Exception:
                      pass
  
                  # identify if this distribution is better
                  if best_sse &gt; sse &gt; 0:
                      best_distribution = distribution
                      best_params = params
                      best_sse = sse
  
          except Exception:
              pass
  
      return (best_distribution.name, best_params)
  
  def make_pdf(dist, params, size=10000):
      """Generate distributions's Probability Distribution Function """
  
      # Separate parts of parameters
      arg = params[:-2]
      loc = params[-2]
      scale = params[-1]
  
      # Get sane start and end points of distribution
      start = dist.ppf(0.01, *arg, loc=loc, scale=scale) if arg else dist.ppf(0.01, loc=loc, scale=scale)
      end = dist.ppf(0.99, *arg, loc=loc, scale=scale) if arg else dist.ppf(0.99, loc=loc, scale=scale)
  
      # Build PDF and turn into pandas Series
      x = np.linspace(start, end, size)
      y = dist.pdf(x, loc=loc, scale=scale, *arg)
      pdf = pd.Series(y, x)
  
      return pdf
  
  # Load data from statsmodels datasets
  data = pd.Series(sm.datasets.elnino.load_pandas().data.set_index('YEAR').values.ravel())
  
  # Plot for comparison
  plt.figure(figsize=(12,8))
  ax = data.plot(kind='hist', bins=50, normed=True, alpha=0.5, color=plt.rcParams['axes.color_cycle'][1])
  # Save plot limits
  dataYLim = ax.get_ylim()
  
  # Find best fit distribution
  best_fit_name, best_fit_params = best_fit_distribution(data, 200, ax)
  best_dist = getattr(st, best_fit_name)
  
  # Update plots
  ax.set_ylim(dataYLim)
  ax.set_title(u'El Niño sea temp.\n All Fitted Distributions')
  ax.set_xlabel(u'Temp (°C)')
  ax.set_ylabel('Frequency')
  
  # Make PDF with best params 
  pdf = make_pdf(best_dist, best_fit_params)
  
  # Display
  plt.figure(figsize=(12,8))
  ax = pdf.plot(lw=2, label='PDF', legend=True)
  data.plot(kind='hist', bins=50, normed=True, alpha=0.5, label='Data', legend=True, ax=ax)
  
  param_names = (best_dist.shapes + ', loc, scale').split(', ') if best_dist.shapes else ['loc', 'scale']
  param_str = ', '.join(['{}={:0.2f}'.format(k,v) for k,v in zip(param_names, best_fit_params)])
  dist_str = '{}({})'.format(best_fit_name, param_str)
  
  ax.set_title(u'El Niño sea temp. with best fit distribution \n' + dist_str)
  ax.set_xlabel(u'Temp. (°C)')
  ax.set_ylabel('Frequency'
  </code></pre>
  
  <p>Update:
  After changing to 'import statsmodels.api as sm'
  I now get..</p>
  
  <pre><code>  File "fitdistr.py", line 106, in &lt;module&gt;
      ax = data.plot(kind='hist', bins=50, normed=True, alpha=0.5, color=plt.rcParams['axes.color_cycle'][1])
    File "/Users/a/.pyenv/versions/3.7.4/lib/python3.7/site-packages/matplotlib/__init__.py", line 835, in __getitem__
      return dict.__getitem__(self, key)
  KeyError: 'axes.color_cycle'
  </code></pre>