Repository: stan-dev/cmdstanpy
Pull request #45 (2019-06-09 16:04:55)

Feature/packaging
#### Submission Checklist

- [x] Run unit tests
- [x] Declare copyright holder and open-source license: see below

#### Summary

- Replaced CmdStan install script `make_cmdstan.sh` with Python script `bin/install_cmdstah.py`- (issue https://github.com/stan-dev/cmdstanpy/issues/41)

- Updated `setup.py` to make `install_cmdstan` available as command line script with package install

- Tested package upload and download via TestPyPI - works on my machine (issue https://github.com/stan-dev/cmdstanpy/issues/17)

- Updated docs to use ScikitLearn templates and styling.

- Edited docs and ipython notebook to reflect current API


#### Copyright and Licensing

Please list the copyright holder for the work you are submitting (this will be you or your assignee, such as a university or company):  Columbia University


By submitting this pull request, the copyright holder is agreeing to license the submitted work under the following licenses:
- Code: BSD 3-clause (https://opensource.org/licenses/BSD-3-Clause)



Patch:
diff --git a/.travis.yml b/.travis.yml
index 89001cd..1c42369 100644
--- a/.travis.yml
+++ b/.travis.yml
@@ -8,10 +8,10 @@ python:
 install:
   - pip install -r requirements.txt
   - pip install flake8 pytest
-  - ./make_cmdstan.sh
+  - ./bin/install_cmdstan
 
 after install:
 
 script:
-  - flake8 --extend-ignore=E127,E201,E202,E203,E231,E252,E999,F841 cmdstanpy
+  - flake8 --extend-ignore=E127,E201,E202,E203,E231,E252,E402,E999,F401,F841 cmdstanpy
   - pytest test
diff --git a/bin/install_cmdstan b/bin/install_cmdstan
new file mode 100755
index 0000000..78868fd
--- /dev/null
+++ b/bin/install_cmdstan
@@ -0,0 +1,166 @@
+#!/usr/bin/env python
+"""
+Download and install a CmdStan release from GitHub.
+Optional command line arguments:
+   -v, --version : version, defaults to latest
+   -d, --dir : install directory, defaults to '~/.cmdstanpy
+"""
+import argparse
+import contextlib
+import os
+import os.path
+import re
+import shutil
+import subprocess
+import sys
+import tarfile
+import urllib.request
+
+
+@contextlib.contextmanager
+def pushd(new_dir):
+    previous_dir = os.getcwd()
+    os.chdir(new_dir)
+    yield
+    os.chdir(previous_dir)
+
+
+def usage():
+    print(
+        """Arguments:
+        -v (--version) :CmdStan version
+        -d (--dir) : install directory
+        -h (--help) : this message
+        """)
+
+
+def install_version(cmdstan_version):
+    with pushd(cmdstan_version):
+        print('Building {} binaries'.format(cmdstan_version))
+        cmd = ['make', 'build']
+        proc = subprocess.Popen(
+            cmd, cwd=None, stdout=subprocess.PIPE, stderr=subprocess.PIPE
+        )
+        proc.wait()
+        stdout, stderr = proc.communicate()
+        if proc.returncode:
+            print('Command "make build" failed')
+            if stderr:
+                print(stderr.decode('utf-8').strip())
+            sys.exit(3)
+        print('Test model compilation')
+        cmd = ['make', os.path.join('examples', 'bernoulli', 'bernoulli')]
+        proc = subprocess.Popen(
+            cmd, cwd=None, stdout=subprocess.PIPE, stderr=subprocess.PIPE
+        )
+        proc.wait()
+        stdout, stderr = proc.communicate()
+        if proc.returncode:
+            print('Failed to compile example model bernoulli.stan')
+            if stderr:
+                print(stderr.decode('utf-8').strip())
+            sys.exit(3)
+    print('Installed {}'.format(cmdstan_version))
+
+
+def is_installed(cmdstan_version):
+    if not os.path.exists(os.path.join(cmdstan_version, 'bin')):
+        return False
+    if os.path.exists(os.path.join(
+            cmdstan_version, 'examples', 'bernoulli', 'bernoulli')):
+        return True
+    if os.path.exists(os.path.join(
+            cmdstan_version, 'examples', 'bernoulli', 'bernoulli.exe')):
+        return True
+    return False
+
+
+def latest_version():
+    try:
+        file_tmp, _ = urllib.request.urlretrieve(
+            'https://api.github.com/repos/stan-dev/cmdstan/releases/latest')
+    except urllib.error.URLError as err:
+        print('Cannot connect to github.')
+        print(err)
+        sys.exit(3)
+    with open(file_tmp, 'r') as fd:
+        response = fd.read()
+        start_idx = response.find('\"tag_name\":\"v') + len('"tag_name":"v')
+        end_idx = response.find('\"', start_idx)
+    return response[start_idx:end_idx]
+
+
+def retrieve_latest_version(version):
+    print('Downloading CmdStan version {}'.format(version))
+    url = 'https://github.com/stan-dev/cmdstan/releases/download/'\
+      'v{0}/cmdstan-{0}.tar.gz'.format(version)
+    try:
+        file_tmp, _ = urllib.request.urlretrieve(url, filename=None)
+    except urllib.error.URLError as err:
+        print('Failed to download CmdStan version {} from github.com'.format(
+            version))
+        print(err)
+        sys.exit(3)
+    print('Download successful, file: {}'.format(file_tmp))
+    try:
+        tar = tarfile.open(file_tmp)
+        tar.extractall(os.getcwd())
+    except Exception as err:
+        print('Failed to unpack download')
+        print(err)
+        sys.exit(3)
+    finally:
+        tar.close()
+    print('Unpacked download as cmdstan-{}'.format(version))
+
+
+def validate_dir(install_dir):
+    if not os.path.exists(install_dir):
+        try:
+            os.makedirs(install_dir)
+        except OSError as e:
+            raise ValueError(
+                'Cannot create directory: {}'.format(install_dir)
+                ) from e
+    else:
+        if not os.path.isdir(install_dir):
+            raise ValueError(
+                'File {} is not a directory: {}'.format(install_dir))
+        try:
+            with open('tmp_test_w', 'w') as fd:
+                pass
+            os.remove('tmp_test_w')  # cleanup
+        except OSError as e:
+            raise ValueError('Cannot write files to directory {}'.format(
+                install_dir)) from e
+
+
+def main():
+    parser = argparse.ArgumentParser()
+    parser.add_argument('--version', '-v')
+    parser.add_argument('--dir', '-d')
+    args = parser.parse_args(sys.argv[1:])
+
+    version = vars(args)['version']
+    if version is None:
+        version = latest_version()
+    print('CmdStan version: {}'.format(version))
+
+    install_dir = vars(args)['dir']
+    if install_dir is None:
+        install_dir = os.path.expanduser(os.path.join('~', '.cmdstanpy'))
+    validate_dir(install_dir)
+    print('Install directory: {}'.format(install_dir))
+
+    cmdstan_version = 'cmdstan-{}'.format(version)
+    with pushd(install_dir):
+        if not os.path.exists(cmdstan_version):
+            retrieve_latest_version(version)
+        if is_installed(cmdstan_version):
+            print('CmdStan version {} already installed'.format(version))
+            sys.exit()
+        install_version(cmdstan_version)
+
+
+if __name__ == "__main__":
+    main()
diff --git a/cmdstanpy/__init__.py b/cmdstanpy/__init__.py
index fc7fa9e..09bb4bb 100644
--- a/cmdstanpy/__init__.py
+++ b/cmdstanpy/__init__.py
@@ -2,7 +2,6 @@ import atexit
 import shutil
 import tempfile
 
-
 STANSUMMARY_STATS = [
     'Mean',
     'MCSE',
@@ -25,3 +24,9 @@ def cleanup_tmpdir():
 
 
 atexit.register(cleanup_tmpdir)
+
+from .cmds import compile_model, sample, summary, diagnose, get_drawset
+from .cmds import save_csvfiles
+from .utils import set_cmdstan_path, cmdstan_path, jsondump, rdump
+from .lib import Model, RunSet
+from ._version import __version__
diff --git a/cmdstanpy/_version.py b/cmdstanpy/_version.py
index b603563..0404d81 100644
--- a/cmdstanpy/_version.py
+++ b/cmdstanpy/_version.py
@@ -1,1 +1,1 @@
-__version__ = "alpha.0.3"
+__version__ = '0.3.0'
diff --git a/cmdstanpy/cmds.py b/cmdstanpy/cmds.py
index d9b16d3..658d78b 100644
--- a/cmdstanpy/cmds.py
+++ b/cmdstanpy/cmds.py
@@ -2,7 +2,6 @@
 First class functions.
 """
 import os
-import os.path
 import platform
 import shutil
 import subprocess
@@ -256,6 +255,15 @@ def diagnose(runset: RunSet) -> None:
     Run cmdstan/bin/diagnose over all output csv files.
     Echo diagnose stdout/stderr to console.
 
+    The diagnose utility reads the outputs of all chains
+    and checks for the following potential problems:
+
+    + Transitions that hit the maximum treedepth
+    + Divergent transitions
+    + Low E-BFMI values (sampler transitions HMC potential energy)
+    + Low effective sample sizes
+    + High R-hat values
+
     :param runset: record of completed run of NUTS sampler
     """
     cmd_path = os.path.join(cmdstan_path(), 'bin', 'diagnose')
@@ -348,10 +356,10 @@ def do_sample(runset: RunSet, idx: int) -> None:
     print('finish chain {}.  '.format(idx + 1))
     with open(transcript_file, 'w+') as transcript:
         if stdout:
-            transcript.write(stdout.decode('ascii'))
+            transcript.write(stdout.decode('utf-8'))
         if stderr:
             transcript.write('ERROR')
-            transcript.write(stderr.decode('ascii'))
+            transcript.write(stderr.decode('utf-8'))
     runset.set_retcode(idx, proc.returncode)
 
 
@@ -367,8 +375,8 @@ def do_command(cmd: str, cwd: str = None) -> str:
     stdout, stderr = proc.communicate()
     if proc.returncode:
         if stderr:
-            msg = 'ERROR\n {} '.format(stderr.decode('ascii').strip())
+            msg = 'ERROR\n {} '.format(stderr.decode('utf-8').strip())
         raise Exception(msg)
     if stdout:
-        return stdout.decode('ascii').strip()
+        return stdout.decode('utf-8').strip()
     return None
diff --git a/cmdstanpy/lib.py b/cmdstanpy/lib.py
index 7a2287c..3f14358 100644
--- a/cmdstanpy/lib.py
+++ b/cmdstanpy/lib.py
@@ -1,5 +1,4 @@
 import os
-import os.path
 import re
 import tempfile
 from typing import Dict, List
diff --git a/cmdstanpy/utils.py b/cmdstanpy/utils.py
index 41d6b16..29f006e 100644
--- a/cmdstanpy/utils.py
+++ b/cmdstanpy/utils.py
@@ -2,24 +2,40 @@
 Utility functions
 """
 import os
-import os.path
 import json
 import numpy as np
 from typing import Dict, TextIO
 
 
+def get_latest_cmdstan(dot_dir: str) -> str:
+    """
+    Given a valid directory path, find all installed CmdStan versions
+    and return highest (i.e., latest) version number.
+    Assumes directory populated via script ``bin/install_cmdstan``.
+    """
+    versions = [
+        name.split('-')[1] for name in os.listdir(dot_dir)
+        if os.path.isdir(os.path.join(dot_dir, name))
+        and name.startswith('cmdstan-')]
+    versions.sort(key=lambda s: list(map(int, s.split('.'))))
+    if len(versions) == 0:
+        return None
+    latest = 'cmdstan-{}'.format(versions[len(versions) - 1])
+    return latest
+
+
 def validate_cmdstan_path(path: str) -> None:
     """
     Validate that CmdStan directory exists and binaries have been built.
     Throws exception if specified path is invalid.
     """
     if not os.path.isdir(path):
-        raise ValueError('no such CmdStan directory {}'.format(path))
+        raise ValueError('no such CmdStan directory {}'.format(
+            path))
     if not os.path.exists(os.path.join(path, 'bin', 'stanc')):
         raise ValueError(
             "no CmdStan binaries found, "
-            "do 'make build' in directory {}".format(path)
-        )
+            "run command line script 'install_cmdstan'")
 
 
 def set_cmdstan_path(path: str) -> None:
@@ -38,8 +54,17 @@ def cmdstan_path() -> str:
     if 'CMDSTAN' in os.environ:
         cmdstan_path = os.environ['CMDSTAN']
     else:
-        cmdstan_path = os.path.expanduser(
-            os.path.join('~', '.cmdstanpy', 'cmdstan'))
+        cmdstan_dir = os.path.expanduser(os.path.join('~', '.cmdstanpy'))
+        if not os.path.exists(cmdstan_dir):
+            raise ValueError(
+                "no CmdStan installation found, "
+                "run command line script 'install_cmdstan'")
+        latest_cmdstan = get_latest_cmdstan(cmdstan_dir)
+        if latest_cmdstan is None:
+            raise ValueError(
+                "no CmdStan installation found, "
+                "run command line script 'install_cmdstan'")
+        cmdstan_path = os.path.join(cmdstan_dir, latest_cmdstan)
     validate_cmdstan_path(cmdstan_path)
     return cmdstan_path
 
diff --git a/cmdstanpy_tutorial.ipynb b/cmdstanpy_tutorial.ipynb
index 2a4709d..3245d92 100644
--- a/cmdstanpy_tutorial.ipynb
+++ b/cmdstanpy_tutorial.ipynb
@@ -12,41 +12,41 @@
     "Given a dataset and a model specification written as a Stan program, the CmdStanPy workflow is:\n",
     "\n",
     "#### Assemble input data as either:\n",
-    "  + a Python `dict` object consisting of key-value pairs where the key corresponds\n",
+    "  + A Python `dict` object consisting of key-value pairs where the key corresponds\n",
     " to Stan data variables and the value is of the correct type and shape.\n",
-    "  + an existing data file on disk in either JSON or Rdump format.\n",
-    "  + CmdStanPy provides methods for converting former to latter.\n",
+    "  + An existing data file on disk in either JSON or Rdump format.\n",
     "\n",
     "#### Compile Stan program: function `compile_model`\n",
-    "  + takes as input a Stan program and produces the corresponding c++ executable.\n",
-    "  + both Stan program and c++ executable exist as on-disk files.\n",
-    "  + returns `Model` object which wraps program, c++ executable filenames or raises exception.\n",
-    "\n",
-    "#### Assemble input data as either:\n",
-    "  + a Python `dict` object consisting of key-value pairs where the key corresponds\n",
-    " to Stan data variables and the value is of the correct type and shape.\n",
-    "  + an existing data file on disk in either JSON or Rdump format.\n",
+    "  + Takes as input a Stan program and produces the corresponding c++ executable.\n",
+    "  + Returns `Model` object which wraps program, c++ executable filenames or raises exception.\n",
+    "  + Both Stan program and c++ executable stored on disk.\n",
     "\n",
     "#### Run sampler: function `sample`\n",
-    "  + invokes Stan's NUTS-HMC sampler to condition model on input data, produces output csv files.\n",
-    "  + creates lightweight object `RunSet` which wraps sampler arguments, per-chain outcome, output files.\n",
-    "  + runs any number of chains - default is 4 chains.\n",
-    "  + validates that all chains run to completion, produce expected number of draws.\n",
-    "  + returns `RunSet` if sample is valid, else raises exception.\n",
+    "  + Invokes Stan's NUTS-HMC sampler to condition model on input data, produces output csv files.\n",
+    "  + Creates lightweight object `RunSet` which wraps sampler arguments, per-chain outcome, output files.\n",
+    "  + Runs any number of chains - default is 4 chains.\n",
+    "  + Validates that all chains run to completion, produce expected number of draws.\n",
+    "  + Returns `RunSet` if sample is valid, else raises exception.\n",
     "\n",
     "#### Assemble fit in-memory: function `get_drawset`\n",
     "  + Sampler output consists of csv files\n",
     "  + `RunSet` method `assemble_sample` reads sampler output csv files into 3-D numpy ndarray (draws, chains, columns)\n",
-    "  + function `get_drawset` flattens chains, returns a pandas DataFrame one row per draw, one column per parameters.\n",
-    "  + argument `params` filters allows list of parameter names.\n",
+    "  + Function `get_drawset` flattens chains, returns a pandas DataFrame one row per draw, one column per parameters.\n",
+    "  + `get_drawset` argument `params` filters DataFrame columns by list of parameter names.\n",
     "\n",
     "#### Summarize posterior fit:  function `summary`\n",
-    "  + invokes CmdStan's `bin/stansummary` utility\n",
-    "  + returns a Pandas DataFrame with estimates of posterior means, standard deviations, Monte-Carlo standard error, effective sample size, and convergence diagnostic statistic for all parameters in the model.\n",
+    "  + Invokes CmdStan's `bin/stansummary` utility\n",
+    "  + Returns a Pandas DataFrame with estimates of posterior means, standard deviations, Monte-Carlo standard error, effective sample size, and convergence diagnostic statistic for all parameters in the model.\n",
     "\n",
     "#### Check for problems: function `diagnose`\n",
-    "  + invokes CmdStan's `bin/diagnose` utility\n",
-    "  + prints output to console\n"
+    "  + Invokes CmdStan's `bin/diagnose` utility\n",
+    "  + Prints output to console\n",
+    "  + Checks for following problems:\n",
+    "    + transitions that hit the maximum treedepth\n",
+    "    + divergent transitions\n",
+    "    + low E-BFMI values (sampler transitions HMC potential energy)\n",
+    "    + low effective sample sizes\n",
+    "    + high R-hat values\n"
    ]
   },
   {
@@ -55,75 +55,92 @@
    "source": [
     "### Installation\n",
     "\n",
-    "* Clone or download CmdStanPy from GitHub repository https://github.com/stan-dev/cmdstanpy\n",
+    "* Install Python package from PyPI or directly from GitHub:\n",
+    "\n",
+    "  + `pip install --upgrade cmdstanpy`\n",
+    "  + `pip install -e git+https://github.com/stan-dev/cmdstanpy`\n",
+    "\n",
     "\n",
     "* CmdStanPy uses CmdStan directly to compile and run Stan programs, therefore CmdStan must be installed locally.\n",
     "\n",
     "  + if you have a working installation of CmdStan, set environment variable `CMDSTAN` to the full path to the top-level CmdStan directory.\n",
     "  \n",
-    "  + if you don't already have CmdStan installed, run script `make_cmdstan.sh` which downloads and compiles the latest release from https://github.com/stan-dev/cmdstan/releases.  By default this installs the latest version of CmdStan in the location `~/.cmdstanpy`.  Flags -d and -v are used to specify the directory and version, respectively.\n"
+    "  + if you don't already have CmdStan installed, run Python script `install_cmdstan` which downloads and compiles the latest release from https://github.com/stan-dev/cmdstan/releases.  By default this installs the latest version of CmdStan in the location `~/.cmdstanpy`.  Flags -d and -v are used to specify the directory and version, respectively.\n"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "### Example 1:  example model `bernoulli.stan`, data `bernoulli.data.json`\n",
+    "### Example 1:  Compile and run CmdStan example model `bernoulli.stan`, data `bernoulli.data.json`\n",
     "\n",
-    "Import relevant classes and methods\n"
+    "Every CmdStan release has an `examples/bernoulli` directory which contains the Stan model and test data files.  In this example we compile the model and run the sampler on the model and data."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Import relevant classes and methods"
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 2,
+   "execution_count": 1,
    "metadata": {},
    "outputs": [],
    "source": [
     "import os\n",
     "import os.path\n",
-    "from cmdstanpy.cmds import compile_model, sample, get_drawset, summary, diagnose"
+    "from cmdstanpy import cmdstan_path, compile_model, sample, get_drawset, summary, diagnose"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Compile model, specify data or data file"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The CmdStanPy directory `test/data` contains models"
+    "The CmdStanPy directory `examples/bernoulli` contains the model and data files."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 3,
+   "execution_count": 5,
    "metadata": {},
    "outputs": [],
    "source": [
-    "datafiles_path = os.path.expanduser(\n",
-    "    os.path.join(\"~\", \"github\", \"stan-dev\", \"cmdstanpy\", \"test\", \"data\"))\n",
-    "stan = os.path.join(datafiles_path, \"bernoulli.stan\")"
+    "bernoulli_path = os.path.join(cmdstan_path(), 'examples', 'bernoulli')\n",
+    "bernoulli_stan = os.path.join(bernoulli_path, 'bernoulli.stan')"
    ]
   },
   {
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Call `compile_model` with path to Stan program, returns `Model` object `bernoulli_model`"
+    "The `compile_model` functions compiles the Stan program and returns a `Model` object."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 4,
+   "execution_count": 7,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "Model(name=bernoulli,  stan_file=\"/Users/mitzi/github/stan-dev/cmdstanpy/test/data/bernoulli.stan\", exe_file=\"/Users/mitzi/github/stan-dev/cmdstanpy/test/data/bernoulli\")\n"
+      "Model(name=bernoulli,  stan_file=\"/Users/mitzi/.cmdstanpy/cmdstan-2.19.1/examples/bernoulli/bernoulli.stan\", exe_file=\"/Users/mitzi/.cmdstanpy/cmdstan-2.19.1/examples/bernoulli/bernoulli\")\n"
      ]
     }
    ],
    "source": [
-    "bernoulli_model = compile_model(stan)\n",
+    "bernoulli_model = compile_model(bernoulli_stan)\n",
     "print(bernoulli_model)"
    ]
   },
@@ -131,12 +148,28 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Input data is either a Python `Dict` with entries corresponding to input data values, or it can be a file in JSON or Rdump format.  If a `Dict` is specified, CmdStanPy writes it to a temp file in JSON format."
+    "Input data is either a Python `Dict` with entries corresponding to input data values, or it can be a file in JSON or Rdump format."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 5,
+   "execution_count": 8,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "bern_json = os.path.join(bernoulli_path, 'bernoulli.data.json')"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "If a `Dict` is specified, CmdStanPy writes it to a temp file in JSON format."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 9,
    "metadata": {},
    "outputs": [],
    "source": [
@@ -147,32 +180,80 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Call `sample` to run NUTS-HMC sampler, returns `PosteriorSample` object `a_sample`"
+    "##### Run the HMC-NUTS sampler on the model and data"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The `sample` function runs the NUTS-HMC sampler and returns a `RunSet` object."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 6,
+   "execution_count": 28,
    "metadata": {},
    "outputs": [
     {
      "name": "stdout",
      "output_type": "stream",
      "text": [
-      "input data tempfile: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmpx89a4v39/tmpbwdt1h2b.json\n",
-      "start chain 1.  start chain 2.  \n",
-      "\n",
+      "input data tempfile: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmpggqfg6d4/tmpjvatu0sj.json\n",
+      "start chain 1.  \n",
+      "finish chain 1.  \n",
+      "start chain 2.  \n",
       "finish chain 2.  \n",
       "start chain 3.  \n",
-      "finish chain 1.  \n",
-      "start chain 4.  \n",
       "finish chain 3.  \n",
+      "start chain 4.  \n",
       "finish chain 4.  \n"
      ]
     }
    ],
    "source": [
-    "bern_fit = sample(bernoulli_model, chains=4, cores=2, data=bern_data)"
+    "bern_fit = sample(bernoulli_model, data=bern_data)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By default, the sample function runs 4 sampler chains.  The `chains` argument specifies the number of chains to run.  The `cores` argument specifies the number of processes to run in parallel."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": 29,
+   "metadata": {},
+   "outputs": [
+    {
+     "name": "stdout",
+     "output_type": "stream",
+     "text": [
+      "input data tempfile: /var/folders/sc/0f0wdc_11_xgjs2v52g20fvr0000gn/T/tmpggqfg6d4/tmp0f38xgj5.json\n",
+      "start chain 1.  start chain 2.  start chain 3.  \n",
+      "\n",
+      "\n",
+      "finish chain 3.  \n",
+      "start chain 4.  \n",
+      "finish chain 1.  \n",
+      "finish chain 2.  start chain 5.  \n",
+      "\n",
+      "finish chain 4.  \n",
+      "finish chain 5.  \n"
+     ]
+    }
+   ],
+   "source": [
+    "bern_fit = sample(bernoulli_model, chains=5, cores=3, data=bern_data)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "##### Summarize or save the results"
    ]
   },
   {
@@ -184,7 +265,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 7,
+   "execution_count": 30,
    "metadata": {},
    "outputs": [
     {
@@ -234,27 +315,27 @@
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>lp__</th>\n",
-       "      <td>-7.292510</td>\n",
-       "      <td>0.023361</td>\n",
-       "      <td>0.816256</td>\n",
-       "      <td>-8.969520</td>\n",
-       "      <td>-6.982420</td>\n",
-       "      <td>-6.749550</td>\n",
-       "      <td>1220.89</td>\n",
-       "      <td>5714.18</td>\n",
-       "      <td>1.00173</td>\n",
+       "      <td>-7.261820</td>\n",
+       "      <td>0.018652</td>\n",
+       "      <td>0.757622</td>\n",
+       "      <td>-8.702790</td>\n",
+       "      <td>-6.978630</td>\n",
+       "      <td>-6.749400</td>\n",
+       "      <td>1649.96</td>\n",
+       "      <td>9970.79</td>\n",
+       "      <td>1.00165</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>theta</th>\n",
-       "      <td>0.255235</td>\n",
-       "      <td>0.003671</td>\n",
-       "      <td>0.122319</td>\n",
-       "      <td>0.080209</td>\n",
-       "      <td>0.241555</td>\n",
-       "      <td>0.473612</td>\n",
-       "      <td>1110.22</td>\n",
-       "      <td>5196.19</td>\n",
-       "      <td>1.00434</td>\n",
+       "      <td>0.250187</td>\n",
+       "      <td>0.002734</td>\n",
+       "      <td>0.119040</td>\n",
+       "      <td>0.082719</td>\n",
+       "      <td>0.237129</td>\n",
+       "      <td>0.467633</td>\n",
+       "      <td>1895.30</td>\n",
+       "      <td>11453.40</td>\n",
+       "      <td>1.00041</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -263,16 +344,16 @@
       "text/plain": [
        "           Mean      MCSE    StdDev        5%       50%       95%    N_Eff  \\\n",
        "name                                                                         \n",
-       "lp__  -7.292510  0.023361  0.816256 -8.969520 -6.982420 -6.749550  1220.89   \n",
-       "theta  0.255235  0.003671  0.122319  0.080209  0.241555  0.473612  1110.22   \n",
+       "lp__  -7.261820  0.018652  0.757622 -8.702790 -6.978630 -6.749400  1649.96   \n",
+       "theta  0.250187  0.002734  0.119040  0.082719  0.237129  0.467633  1895.30   \n",
        "\n",
-       "       N_Eff/s    R_hat  \n",
-       "name                     \n",
-       "lp__   5714.18  1.00173  \n",
-       "theta  5196.19  1.00434  "
+       "        N_Eff/s    R_hat  \n",
+       "name                      \n",
+       "lp__    9970.79  1.00165  \n",
+       "theta  11453.40  1.00041  "
       ]
      },
-     "execution_count": 7,
+     "execution_count": 30,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -290,7 +371,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 8,
+   "execution_count": 31,
    "metadata": {},
    "outputs": [
     {
@@ -309,63 +390,46 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `get_drawset` function returns a pandas.DataFrame, one draw per row.  It contains all columns from the sampler's csv output file, i.e., it contains both the sampler state and the values for all parameter, transformed parameter, and generated quantities variables."
+    "The `get_drawset` function returns a pandas.DataFrame, one draw per row."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 9,
+   "execution_count": 32,
    "metadata": {},
-   "outputs": [
-    {
-     "data": {
-      "text/plain": [
-       "((4000, 8),\n",
-       " Index(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__',\n",
-       "        'divergent__', 'energy__', 'theta'],\n",
-       "       dtype='object'))"
-      ]
-     },
-     "execution_count": 9,
-     "metadata": {},
-     "output_type": "execute_result"
-    }
-   ],
+   "outputs": [],
    "source": [
-    "bern_drawset = get_drawset(bern_fit)\n",
-    "bern_drawset.shape, bern_drawset.columns"
+    "bern_drawset = get_drawset(bern_fit)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "By default, `get_drawset` returns a DataFrame which contains all columns from the sampler's csv output file, i.e., it contains both the sampler state and the values for all parameter, transformed parameter, and generated quantities variables."
    ]
   },
   {
    "cell_type": "code",
-   "execution_count": 11,
+   "execution_count": 33,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "<matplotlib.axes._subplots.AxesSubplot at 0x114636d68>"
+       "((5000, 8),\n",
+       " Index(['lp__', 'accept_stat__', 'stepsize__', 'treedepth__', 'n_leapfrog__',\n",
+       "        'divergent__', 'energy__', 'theta'],\n",
+       "       dtype='object'))"
       ]
      },
-     "execution_count": 11,
+     "execution_count": 33,
      "metadata": {},
      "output_type": "execute_result"
-    },
-    {
-     "data": {
-      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xd4XPWd7/H3V11Wsy3JTXK3MDG2wVhUQwIkEEM2wAZCyN4UshCWQMomu7mXze6Tm83e++wme5M8yYZNQkghnRRCHAIpEIdisI1t3LvlJjdJtmwV26rf+8eMhCyrjMqZMyN9Xs8zD2fO/GbmgzzSd875/X7nZ+6OiIgIQErYAUREJHGoKIiISCcVBRER6aSiICIinVQURESkk4qCiIh0UlEQEZFOgRUFM8sys9VmtsHMtpjZv/bQ5h4zqzaz9dHbfUHlERGR/qUF+NpNwA3u3mBm6cDLZvasu6/s1u4Jd/9ogDlERCRGgRUFj0yVbojeTY/ehjx9uqioyGfMmDHUlxERGVXWrl1b4+7F/bUL8kgBM0sF1gJzgEfcfVUPze4wszcDO4FPuvvBvl5zxowZrFmzZvjDioiMYGa2P5Z2gXY0u3ubu18ClAKXm9n8bk1+C8xw94XAn4DHe3odM7vfzNaY2Zrq6uogI4uIjGpxGX3k7ieB5cDSbvuPu3tT9O5jwOJenv+ou5e7e3lxcb9HPyIiMkhBjj4qNrOx0e1s4EZge7c2k7vcvRXYFlQeERHpX5B9CpOBx6P9CinAz939aTP7PLDG3ZcBHzezW4FW4ARwT4B5RESkH5Zs6ymUl5e7OppFRAbGzNa6e3l/7TSjWUREOqkoiIhIp0DnKYjEQ2tbO7/bdISjp85yy4LJTB0/JuxIIklLRUGS2pnmNu753mpW7T0BwFee28n37rmcq2YXhpxMJDnp9JEkLXfnUz9fz2v7TvDFOxfy0v+8nqnjxvCxn66jtrE57HgiSUlFQZLWsg2HeXbzUT799gu5q3wqU8eP4at3L+J4YzOPvlQRdjyRpKSiIEmpqbWNLzy7nYWlBdz/5lmd++dNyeeW+ZP50av7OXWmJcSEIslJRUGS0s9fO8jhU2f59Nvnkppi5zx2/5tnUd/UytMbD4eUTiR5qShI0jnb0sbXl+/mshnjuGZO0XmPLywtYM6EXH7zuoqCyECpKEjS+dnqAxyra+KTN16AmZ33uJlx+yVTWL3vBIdPngkhoUjyUlGQpHK2pY1H/rKHK2eN5+rZ5x8ldFg6fxIAy3dUxSuayIigoiBJ5Ucr91Nd38Qn33ZBn+1mF+dSOi6b5du1/obIQKgoSNI43dzKN1/YwzVzirhiVt+T08yM6+YWs2J3DU2tbXFKKJL8VBQkaXz/lX3UNDTzyRvLYmp/bVkxZ1ra2Fh5KuBkIiOHioIkhVOnW/jmX/Zww4UTWDx9fEzPuWxGpN3q6CUwRKR/KgqSFL714h7qzrbyjzfNjfk543MyKJuQq6IgMgAqCpLwqurP8r0V+7j14inMm5I/oOdePnM86/bX0taeXItJiYRFRUES3tee30VLWzufurHvEUc9uXzmeOqbWtl2pC6AZCIjj4qCJLQ91Q38dPVB/uaKacwoyhnw8xdPHwfA6wdqhzuayIikoiAJ7Yu/305WWgoff2tsI466KxmbTWFOBhs0AkkkJoEVBTPLMrPVZrbBzLaY2b/20CbTzJ4ws91mtsrMZgSVR5LPmn0n+MOWYzzwltkU5WYO6jXMjIWlBWysPDnM6URGpiCPFJqAG9z9YuASYKmZXdmtzb1ArbvPAb4CfCHAPJJE3J3/eHY7E/IyuffamUN6rYWlY9ld1UBjU+swpRMZuQIrCh7REL2bHr11HwJyG/B4dPuXwFutpyucyaizau8J1uyv5WM3zGFMxtBWjb14agHtDpsP6RSSSH8C7VMws1QzWw9UAX9y91XdmpQABwHcvRU4BWhxXeEbf9lDUW4G7y6fOuTXWlg6FkAzm0ViEGhRcPc2d78EKAUuN7P5g3kdM7vfzNaY2Zrqal3gbKTbcvgUL+ys5kNLZpKVnjrk1yvKzaRkbDYb1K8g0q+4jD5y95PAcmBpt4cOAVMBzCwNKACO9/D8R9293N3Li4uLg44rIfvRygNkp6fyviunD9trLiwtYJNOH4n0K8jRR8VmNja6nQ3cCGzv1mwZ8MHo9p3An91dU09HsdPNrfx2w2FuWTCZguz0YXvd+SUF7D9+Wus2i/QjyCOFycByM9sIvEakT+FpM/u8md0abfMdoNDMdgOfAh4OMI8kgWc3HaWhqZW7ykuH9XXnlxQAkVNTItK7oQ3r6IO7bwQW9bD/s122zwLvDiqDJJ+n1h9i2vgxXD4ztiuhxmp+9JpJWw7V9blim8hopxnNkjBOnWnh1T3HuXn+pB7XXh6KwtxMphRkqV9BpB8qCpIwlm+vorXduemiSYG8/vySAjbr9JFIn1QUJGH8YctRJuRlsmjq2EBef35JAXtrGmnQzGaRXqkoSEJobWvn5V01vPVNE0hJCWZS+4KSAtxh62FdRlukNyoKkhA2HjpFfVMr18wJbh7KRSWRzmb1K4j0TkVBEsKKXTUAXDU7uKucTMjLYmJ+JltUFER6paIgCWHFnhoumpLP+JyMQN9n/hTNbBbpi4qChO5sSxvr9p9kyZzg5w/MLylgT3UDp5vV2SzSExUFCd3GylM0t7Vz+YzhnbDWk/klkctoa81mkZ6pKEjoOtZPXjQtmKGoXS2IXu5i8yEVBZGeqChI6NYdqGV64RgKB7nk5kBMzM+kKDdD/QoivVBRkFC5O+sOnOTSaePi8n5mFpnZrKIg0iMVBQlVZe0ZquubuDQOp446zJ9SwK6qBs62tMXtPUWShYqChGpdZ39CfI4UINLZ3NbubD9aH7f3FEkWKgoSqvUHT5KdnsqFk/Li9p7zNbNZpFcqChKqLYfqmDcln7TU+H0US8ZmM25MumY2i/RARUFC097ubD1Sx7zJ+XF9347OZh0piJxPRUFCc7D2NA1NrVw0Jb5FASL9CjuP1dPUqs5mka5UFCQ0W6KXsJ4XRlGYUkBLm7PzaEPc31skkakoSGi2Hq4jNcW4YGL8Opk7dM5s1kpsIudQUZDQbDl8irIJuWSlp8b9vaeOzyY/K039CiLdBFYUzGyqmS03s61mtsXMPtFDm+vM7JSZrY/ePhtUHkk8Ww7Hv5O5Q0dns0YgiZwrLcDXbgX+wd3XmVkesNbM/uTuW7u1e8nd/yrAHJKAahqaqKpvCqU/ocP8kgK+/8o+WtraSY/jkFiRRBbYb4K7H3H3ddHtemAbUBLU+0ly2XksMpt4bhwnrXU3v6SA5tZ2dh1TZ7NIh7h8PTKzGcAiYFUPD19lZhvM7FkzuygeeSR8HX+Iw+hk7jA/epSii+OJvCHwomBmucCvgL939+4XsV8HTHf3i4H/Ap7q5TXuN7M1Zramuro62MASFzuP1ZOflcaEvOAvl92bGYU55GamsfHQydAyiCSaQIuCmaUTKQg/dvcnuz/u7nXu3hDdfgZIN7Pz1mR090fdvdzdy4uLi4OMLHGyq6qBsol5mFloGVJSjEXTxrJmX21oGUQSTZCjjwz4DrDN3b/cS5tJ0XaY2eXRPMeDyiSJwd3ZdayeCybmhh2FK2aOZ/vRek6ebg47ikhCCHL00RLg/cAmM1sf3fcZYBqAu38TuBP4iJm1AmeAu93dA8wkCaCmoZna0y2UTQivP6HD5TMLAVi99wQ3XTQp5DQi4QusKLj7y0Cf5wbc/evA14PKIIlpV3TkUZidzB0unlpARlqKioJIlAZnS9ztqoqMPCpLgNNHmWmpLJo6llV7T4QdRSQhqChI3CXCyKOurpg5ni2HT1F/tiXsKCKhU1GQuNt1rIELQh551NUVswppd1i7X6OQRFQUJK7cnZ1V9ZQlQH9Ch0XTxpKearyyRwPfRFQUJK5qGpo5ebolIYajdhiTkcbi6eN4aVdN2FFEQqeiIHG1qyoy8mjOhMQpCgDXlhWz7Ugd1fVNYUcRCZWKgsRVRXUjALOLE60oRCbSr9itowUZ3VQUJK721jSSnZ7KpPyssKOc46IpBYwbk86Lu3RtLRndVBQkriqqG5hRlENKSmKMPOqQmmIsmVPEy7tq0KR6Gc1UFCSu9tY0Mqs4J+wYPbq2rIiq+iZ2an0FGcVUFCRumlvbOVh7hllFiVkUrimLXIH3JZ1CklFMRUHi5sCJ07S1e8IeKZSMzWZ2cQ4vamiqjGIqChI3FdWR0zKzihJr5FFX15YVs6riOGdb2sKOIhIKFQWJm4qayHDUmQl6pABw/YUTaGpt51XNbpZRSkVB4mZvdSNFuZnkZ6WHHaVXV8wcT3Z6Kst3VIUdRSQUKgoSNxU1DQnbydwhKz2VJXMK+fP2Kg1NlVFJRUHiJpGHo3Z13dwJVNaeYU+1hqbK6KOiIHFx6kwLNQ3NzEzwIwWI9CsALN+uoaky+qgoSFzsjXYyz0qwax71pGRsNnMn5qlfQUYlFQWJi87hqElw+gjguguLeW3fCa3GJqNOYEXBzKaa2XIz22pmW8zsEz20MTP7mpntNrONZnZpUHkkXBXVjaSmGFPHjQk7SkyunzuBljbXVVNl1AnySKEV+Ad3nwdcCTxkZvO6tbkZKIve7ge+EWAeCdHemkamjR9DRlpyHJwunj6OvKw09SvIqBPYb6i7H3H3ddHtemAbUNKt2W3ADzxiJTDWzCYHlUnCs6e6ISk6mTukp6bw5rJilu/Q0FQZXeLytc3MZgCLgFXdHioBDna5X8n5hUOSXHu7s+94Y8LPUejuurnFVNU3seVwXdhRROIm8KJgZrnAr4C/d/dB/XaZ2f1mtsbM1lRX63A+2RypO8vZlvakGHnU1XVzI0NT/7xdo5Bk9Ai0KJhZOpGC8GN3f7KHJoeAqV3ul0b3ncPdH3X3cncvLy4uDiasBGZvdAnOZDp9BFCcl8ml08by+81Hw44iEjcxFQUze9LM3mFmMRcRMzPgO8A2d/9yL82WAR+IjkK6Ejjl7kdifQ9JDhU1keGos5NkOGpXtyyYzNYjdeyLzrMQGeli/SP/38DfALvM7D/MbG4Mz1kCvB+4wczWR2+3mNkDZvZAtM0zQAWwG/g28OAA80sSqKhuJCcjleK8zLCjDNjNCyLjHp7ZrO8qMjqkxdLI3Z8DnjOzAuC90e2DRP6Q/8jdz5vh4+4vA30uxOuRYR0PDTi1JJWKmkZmFecSOXhMLiVjs7lk6lie2XSEB6+bE3YckcAN5HRQIXAPcB/wOvBV4FLgT4EkkxGjorohaWYy9+SWBZPYfKiOA8dPhx1FJHCx9in8GngJGAO8091vdfcn3P1jQHINKZG4OtvSxqGTZ5Kuk7mrW6KnkH6z/rwxECIjTqxHCt9293nu/u8dHcFmlgng7uWBpZOkt//4adyT40J4vSkdN4arZhXyi7WVtLdrIpuMbLEWhf/Tw75XhzOIjEx7azrWZU7eIwWA91w2lQMnTrNyr5bplJGtz45mM5tEZIZxtpkt4o2O43wip5JE+rQnSecodLd0/iTyfpPGz187yNWzi8KOIxKY/kYfvZ1I53Ip0HWuQT3wmYAyyQiyt6aRifmZ5GTGNNAtYWWlp3L7JSU8seYg/9LQRFFu8g2vFYlFn6eP3P1xd78euMfdr+9yu7WXGcoi59hb08isouTtT+jqniUzaG5t50cr94cdRSQwfRYFM3tfdHOGmX2q+y0O+STJVVQ3MDOJh6N2Nbs4lxsunMAPX93P2Za2sOOIBKK/juaO3+ZcIK+Hm0ivahubqT3dkvSdzF3dd81Mjjc2s2z94bCjiASizxO97v6t6H//NT5xZCTZe3xkdDJ3ddXsQi6clMdjL1fw7vLSpJylLdKXWCevfdHM8s0s3cyeN7PqLqeWRHpUER15lMxzFLozM+69ZiY7jzXwyh4NT5WRJ9Z5CjdF10L4K2AfMAf4dFChZGTYW9NAWopROi477CjD6p0XT6EwJ4PvrdgbdhSRYRdrUeg4zfQO4BfufiqgPDKCdKzLnJ6aHOsyxyorPZW/uWIaz2+vYv9xXVJbRpZYf1ufNrPtwGLgeTMrBs4GF0tGgorqxqS+EF5f3nfldFLNePwVDU+VkSWmouDuDwNXA+XRy2Q3ArcFGUySW8e6zCOpk7mriflZ3LxgMr9ce1DDU2VEGchx/YXAe8zsA8CdwE3BRJKRoGNd5pkjZOJaT+4qL6XubKvWcJYRJaZrD5jZD4HZwHqg42uRAz8IKJckuWRdl3kgrp5dxKT8LH61trLz8toiyS7WC9KUA/OiK6WJ9CuZ12WOVWqKcfuiEr79UgUnGpsZn5MRdiSRIYv19NFmYFKQQWRkSeZ1mQfirxZOpq3deW7bsbCjiAyLWI8UioCtZrYaaOrY6e63BpJKkt7emkZmFueM+Bm/F03Jp2RsNn/ccpS7yqeGHUdkyGItCp8LMoSMPHtrGrl46tiwYwTOzLhx3kR+svoAjU2tSX+JcJFYh6S+QGQmc3p0+zVgXV/PMbPvmlmVmW3u5fHrzOyUma2P3j47wOySoJpa26isPT2iLoTXl7dfNInm1nZe2lUddhSRIYv12kcfBn4JfCu6qwR4qp+nfR9Y2k+bl9z9kujt87FkkcR34Php2p0RO3Gtu/IZ48jJSOXFXTVhRxEZslg7mh8ClgB1AO6+C5jQ1xPc/UXgxJDSSVKqqBn5w1G7Sk9N4arZRby4sxoN0JNkF2tRaHL35o47ZpZGZJ7CUF1lZhvM7Fkzu6i3RmZ2v5mtMbM11dU6RE90FaNgjkJ3b76giMraM+w/fjrsKCJDEmtReMHMPgNkm9mNwC+A3w7xvdcB0939YuC/6ON0lLs/6u7l7l5eXFw8xLeVoO2taaA4L5O8rPSwo8TNtWWRz6X6FSTZxVoUHgaqgU3A3wHPAP8ylDd29zp3b4huPwOkm1nRUF5TEsPempF7zaPezCgcw9Tx2bywU/0KktxiGj/n7u1m9hTwlLsPy1chM5sEHHN3N7PLiRQorVoyAuytaeRtb5oYdoy4MjOumVPM0xsO09rWTtoIu1y4jB59fnIt4nNmVgPsAHZEV13rd/iomf0UeBWYa2aVZnavmT1gZg9Em9wJbDazDcDXgLt1GY3kd+pMCzUNzaNm5FFXV88upL6plc2H68KOIjJo/R0pfJLIqKPL3H0vgJnNAr5hZp9096/09kR3f29fL+zuXwe+PsC8kuD2VEeueTSSr47am6tmFwKwYncNl4yCiXsyMvV3jPt+4L0dBQHA3SuA9wEfCDKYJKfdVZGiMGfC6CsKRbmZXDgpj1f2qF9Bkld/RSHd3c/7hEf7FUbP0BKJ2Z6qBjJSU5g6wtZljtXVs4tYs69WC+9I0uqvKDQP8jEZpXZXNTCzKGfUdrRePbuQptZ21h2oDTuKyKD095t7sZnV9XCrBxbEI6Akl93VDaPy1FGHK2aNJzXFeHWPBtJJcuqzKLh7qrvn93DLc3edPpJznG1p4+CJ08wexUUhLyudBSUFrNitfgVJTqPzGF8CsbemkXYfnZ3MXS2ZU8iGylM0NLWGHUVkwFQUZNh0jjwqHt1F4erZRbS1O6v36hSSJB8VBRk2u6saMBs9l8zuzeLp48hIS+GV3SoKknxUFGTY7K5uYNr4MWSlp4YdJVRZ6amUTx/HCnU2SxJSUZBhs6eqYdSfOupw9exCth2p40SjRm5LclFRkGHR2tZORXXjqO9k7nD1nMgFfzU0VZKNioIMi4O1Z2huax/Vw1G7WlhSQG5mmi55IUlHRUGGxWi+5lFP0lJTuGLmeF7RkYIkGRUFGRYqCue7anYhe2saOXzyTNhRRGKmoiDDYndVAxPyMskfRUtw9mdJtF9BRwuSTFQUZFiM9mse9WTuxDzG52Twii55IUlERUGGzN0jw1FVFM6RkmJcNbuQV/YcR4sKSrJQUZAhq6w9Q0NTK3Mn5YUdJeFcPbuQo3VnqahpDDuKSExUFGTIdhytB+BCFYXzvLmsGIDl26tCTiISGxUFGbIdxyJFoWyiikJ3U8ePYd7kfJ7dfDTsKCIxCawomNl3zazKzDb38riZ2dfMbLeZbTSzS4PKIsHacbSekrHZGnnUi6XzJ7F2fy1VdWfDjiLSryCPFL4PLO3j8ZuBsujtfuAbAWaRAO04Wq/+hD7cPH8SAH/YeizkJCL9C6wouPuLwIk+mtwG/MAjVgJjzWxyUHkkGM2t7eypblBR6MOcCbnMKs7h95uPhB1FpF9h9imUAAe73K+M7pMksremkdZ2Z676E3plZtw8fxIrK07oqqmS8JKio9nM7jezNWa2prq6Ouw40sX2o3UAOlLoxzsWTKGt3fndJh0tSGILsygcAqZ2uV8a3Xced3/U3cvdvby4uDgu4SQ2O4/Vk5ZizNY6Cn160+Q85k7M46nXe/yIiySMMIvCMuAD0VFIVwKn3F1fo5LMjqP1zCzKISMtKQ46Q2Nm3L6ohLX7azlw/HTYcUR6FeSQ1J8CrwJzzazSzO41swfM7IFok2eACmA38G3gwaCySHC2a+RRzG67ZAoAT63X0YIkrrSgXtjd39vP4w48FNT7S/AamlqprD3D3ZdN7b+xMGVsNlfOGs9Trx/iYzfMwczCjiRyHh3zy6DtjM5kvkAjj2L214tKqKhpZEPlqbCjiPRIRUEG7Y1rHuWHnCR5LJ0/mYy0FHU4S8JSUZBB23G0njEZqZSOyw47StIoyE7nbW+awLINh2lpaw87jsh5VBRk0HYcreeCiXmkpOjc+EC8a1EpJxqbeWGH5txI4lFRkEFxd3Ycq9dM5kF4y9xiCnMy+NW6yrCjiJxHRUEGpbqhiRONzRqOOgjpqSnceskUnt9WxcnTuuyFJBYVBRkULawzNHdcWkpzWzu/3aj5mpJYVBRkULYejlzz6E2TNfJoMC6aks/ciXk8qVNIkmBUFGRQth2pY3JBFuNyMsKOkpTMjHddWsLrB06yp7oh7DginVQUZFC2Hqljno4ShuT2RSWkGPx6neYsSOJQUZABO9vSxp7qRuZNUVEYion5WVxTVsyvXz9Ee7uHHUcEUFGQQdh1rIG2dld/wjC449ISDp08w8q9x8OOIgKoKMggbD0SuW6PTh8N3U3zJpGbmcav1uoUkiQGFQUZsK2H68jJSGXa+DFhR0l62RmpvGPBZJ7dfITTza1hxxFRUZCB23qkjjdNztflLYbJuy4t4XRzG7/ffDTsKCIqCjIw7e3OtiP16k8YRpfNGM/U8dk8qVFIkgBUFGRAKmvP0NDUqpFHwyglxfjrRaWs2FPD4ZNnwo4jo5yKggyIOpmDccelJbhrqU4Jn4qCDMjWw3WkGLoQ3jCbXphD+fRx/GptJZGVakXCoaIgA7L1SD2zinPJSk8NO8qIc8fiUvZUN7JRS3VKiFQUZEC26fIWgbllQWSpTq2zIGEKtCiY2VIz22Fmu83s4R4ev8fMqs1sffR2X5B5ZGhONDZz6OQZdTIHpCA7nZvmTWTZhsM0t2qpTglHYEXBzFKBR4CbgXnAe81sXg9Nn3D3S6K3x4LKI0O36VDktMbC0oKQk4xcd1xaysnTLfx5+7Gwo8goFeSRwuXAbnevcPdm4GfAbQG+nwRs48GTAMwvUVEIyrVlRUwpyOK7K/aFHUVGqSCLQglwsMv9yui+7u4ws41m9kszm9rTC5nZ/Wa2xszWVFdrsfOwbDx0ilnFOeRnpYcdZcRKS03hvmtnsXrvCdbuPxF2HBmFwu5o/i0ww90XAn8CHu+pkbs/6u7l7l5eXFwc14Dyhk2Vp1ioo4TA3X35VMaOSecbf9kTdhQZhYIsCoeArt/8S6P7Orn7cXdvit59DFgcYB4Zgqq6sxytO8uC0rFhRxnxxmSk8aGrZ/LctipeP1AbdhwZZYIsCq8BZWY208wygLuBZV0bmNnkLndvBbYFmEeGQJ3M8XXftTMpzsvk357eqslsEleBFQV3bwU+CvyByB/7n7v7FjP7vJndGm32cTPbYmYbgI8D9wSVR4ZmQ+UpUiyy4LwELyczjU/fNJd1B06ybMPhsOPIKJIW5Iu7+zPAM932fbbL9j8B/xRkBhkemypPUjYhjzEZgX5kpIs7Fpfyo1X7+dyyLVw9u4jivMywI8koEHZHsySB9nZn/cGTOnUUZ6kpxpfefTGNzW3805MbdRpJ4kJFQfpVUdNA7ekWLpsxPuwoo07ZxDz+19ILeW5bFf+t0UgSByoK0q/X9kVGwJTPGBdyktHpb5fM4NaLp/D//riD57ZqprMES0VB+rVmXy2FORnMLMoJO8qoZGZ84Y6FXDQln4//7HUNU5VAqShIv9bsP0H5jHGYaU3msGRnpPLdD15GUW4m93zvNXYcrQ87koxQKgrSp6r6s+w/fpry6epPCNuE/Cx+dO8VZKal8P7vrOLA8dNhR5IRSEVB+vTaXvUnJJJphWP44b1X0NTazt2PvkpFdUPYkWSEUVGQPr20q5q8rDQW6JpHCWPupDx+8uFIYbjrWyvZdqQu7EgygqgoSK/cnRd3VrNkdhFpqfqoJJKLphTwxN9dRWoK3PmNV3h205GwI8kIod906dWe6gYOnzrLmy/QlWkT0ZwJuTz10BLKJubxkR+v45NPrKeq/mzYsSTJqShIr17YWQNEFn6RxDS5IJsn/u5KPnr9HJ7eeJhrv7Cch3+1kc3RCxiKDJQuZCO9+uOWo5RNyGXq+DFhR5E+ZKal8o9vn8sdi0t59MUKnlxXyc9eO8jC0gLuvmwat14yhdxM/apLbHSkID2qqjvL6n0neMfCyf03loQwsyiHf3/XAlZ/5m187p3zaGpp5zO/3sTl//c5PrdsC6fOtIQdUZKAvj5Ij57dfBR3eMcCFYVkUzAmnXuWzOSDV8/g9YMn+fHKA/zg1X38btMR/vPOhVw3d0LYESWB6UhBevSb9Ye4YGIuZRPzwo4ig2RmXDptHF+662KWffQaCnMy+ND3X+O/nt9Fe7uuuCo9U1GQ82w9XMe6Aye5q3xq/40lKcwvKeDXDy7htoun8KU/7eT+H66l7qxOJ8n5VBTkPD9cuY/MtBTuXFwadhSAhlo3AAAJRUlEQVQZRtkZqXzlPZfwv985j+U7qrj9kRXsrtKMaDmXioKco7L2NL9cW8mdi0sZOyYj7DgyzMyMDy2ZyY/vu4JTp1u4/ZEV/GTVAZ1Okk4qCnKOL/9xJ4bx0PVzwo4iAbpyViFPf/wa5pfk85lfb+L2/17BM5uO0NzaHnY0CZlGH0mnF3ZW8+Trh3jo+tlMGZsddhwJ2OSCbH764St5ct0hvvbnXTz443XkZqbxlguKubasiCVzijRHZRQKtCiY2VLgq0Aq8Ji7/0e3xzOBHwCLgePAe9x9X5CZpGcV1Q184mevUzYhl4/dUBZ2HIkTM+OOxaXcvqiEF3dW88etR3l+WxW/i15LaXrhGJbMKeKaOUVcNauQcTk6pTjSBVYUzCwVeAS4EagEXjOzZe6+tUuze4Fad59jZncDXwDeE1Qm6dkru2v42E9fJ8WMxz5YTlZ6atiRJM5SU4zrL5zA9RdOwN3ZU93Ay7tqeHn3cZatP8xPVh3ADOZNzueKmYVcMWs880sKmJyfRUqKFl8aScw9mA4mM7sK+Jy7vz16/58A3P3fu7T5Q7TNq2aWBhwFir2PUOXl5b5mzZpAMo8WZ1vaOHzyDGv317Jsw2Fe2lXDrKIcHvtgObOKc8OOJwmmpa2djZWneHlXDa9W1PD6gZM0RfsestJTmFGYw+SCLIrzMiO33EyK87IYl5NOTkYaYzJSGZOZRnqKkZJipJiRYnRup5phFilMHY9plb/hZ2Zr3b28v3ZBnj4qAQ52uV8JXNFbG3dvNbNTQCFQM9xhXthZzb89HTlI6ag551QeP3+za23q2rZjt3fZ27WM9VTSenrPc54TfcR7yHHu/l7es4/MXdu1tzv1Ta2d9yfmZ/Lpt8/lb5fMJDtDRwhyvvTUFBZPH8fi6eP4BGU0tbaxqfIU24/WU1HdyL7jjVTVn2XL4TqONzbTNgwjmcwg1aJFIoVosXijmEQKiZHa9bGUrs8xRmJZec9lU7nv2lmBvkdSdDSb2f3A/QDTpk0b1GvkZqYxt+vsXDvnPx3v0/1hun5h6antOR+8c9par8+3Htp13W+9vKj1mDnW93xjZ1FuBpMKsllQUsAFE3P1rUwGJDMtlfIZ4ymfcf4SrW3tTu3pZqrrm6g93cyZ5jZON7dxprmNlvZ22j3yxaTdnbZ2xx3aPHI/sp/ofo/u79qeSDt/475HX6fd6bI/8t+RqCg3M/D3CLIoHAK6Toktje7rqU1l9PRRAZEO53O4+6PAoxA5fTSYMB3fdEQkOKkpRlFuZlz+eEkwgpyn8BpQZmYzzSwDuBtY1q3NMuCD0e07gT/31Z8gIiLBCuxIIdpH8FHgD0SGpH7X3beY2eeBNe6+DPgO8EMz2w2cIFI4REQkJIH2Kbj7M8Az3fZ9tsv2WeDdQWYQEZHY6TIXIiLSSUVBREQ6qSiIiEgnFQUREemkoiAiIp0Cu/ZRUMysGtg/wKcVEcClM4ZBouYCZRusRM2WqLlA2QZroNmmu3txf42SrigMhpmtieVCUPGWqLlA2QYrUbMlai5QtsEKKptOH4mISCcVBRER6TRaisKjYQfoRaLmAmUbrETNlqi5QNkGK5Bso6JPQUREYjNajhRERCQGI7IomNl4M/uTme2K/rfXhRTMLN/MKs3s64mQy8wuMbNXzWyLmW00s0DXrDazpWa2w8x2m9nDPTyeaWZPRB9fZWYzgswzwGyfMrOt0Z/T82Y2PRFydWl3h5m5mcVt9Eos2czsrujPbYuZ/SRRspnZNDNbbmavR/9Nb4lTru+aWZWZbe7lcTOzr0VzbzSzS+ORK8Zs/yOaaZOZvWJmFw/5Td19xN2ALwIPR7cfBr7QR9uvAj8Bvp4IuYALgLLo9hTgCDA2oDypwB5gFpABbADmdWvzIPDN6PbdwBNx+jeMJdv1wJjo9kfikS2WXNF2ecCLwEqgPIF+ZmXA68C46P0JCZTtUeAj0e15wL44ZXszcCmwuZfHbwGeJbKQ4ZXAqnjkijHb1V3+LW8ejmwj8kgBuA14PLr9OHB7T43MbDEwEfhjouRy953uviu6fRioAvqdcDJIlwO73b3C3ZuBn0Uz9pb5l8BbLT7rd/abzd2Xu/vp6N2VRFb3Cz1X1L8BXwDOxiHTQLJ9GHjE3WsB3L0qgbI5kB/dLgAOxyOYu79IZD2X3twG/MAjVgJjzWxyImRz91c6/i0Zpt+BkVoUJrr7kej2USJ/+M9hZinAl4B/TKRcXZnZ5US+Ve0JKE8JcLDL/crovh7buHsrcAooDCjPQLN1dS+Rb3NB6zdX9PTCVHf/XRzydBXLz+wC4AIzW2FmK81saQJl+xzwPjOrJLIOy8fiE61fA/0shmVYfgcCXWQnSGb2HDCph4f+uesdd3cz62mI1YPAM+5eOZxffIchV8frTAZ+CHzQ3duHLeAIZGbvA8qBtyRAlhTgy8A9IUfpTRqRU0jXEflW+aKZLXD3k6Gmingv8H13/5KZXUVkVcb5+vz3z8yuJ1IUrhnqayVtUXD3t/X2mJkdM7PJ7n4k+se1p0Pkq4BrzexBIBfIMLMGd++14zBOuTCzfOB3wD9HD1eDcgiY2uV+aXRfT20qzSyNyGH98QAzDSQbZvY2IgX3Le7elAC58oD5wF+iXzYmAcvM7FZ3XxNyNoh8y13l7i3AXjPbSaRIvJYA2e4FlgK4+6tmlkXk+j7xOsXVm5g+i2Exs4XAY8DN7j703814dZjE8wb8J+d26H6xn/b3EJ+O5n5zETld9Dzw93HIkwZUADN5o/Pvom5tHuLcjuafx+nfMJZsi4icWiuL42er31zd2v+F+HU0x/IzWwo8Ht0uInJapDBBsj0L3BPdfhORPgWL089uBr135r6DczuaV8fr8xZDtmnAbuDqYXu/eP7PxfGHWBj9w7oLeA4YH91fDjzWQ/t4FYV+cwHvA1qA9V1ulwSY6RZgZ/SP6z9H930euDW6nQX8IvrBWw3MiuO/Y3/ZngOOdfk5LUuEXN3axq0oxPgzMyKnt7YCm4C7EyjbPGBFtGCsB26KU66fEhnl10LkSOpe4AHggS4/s0eiuTfF+d+zv2yPAbVdfgfWDPU9NaNZREQ6jdTRRyIiMggqCiIi0klFQUREOqkoiIhIJxUFERHppKIgIiKdVBRERKSTioKIiHT6//+biJG6YTPAAAAAAElFTkSuQmCC\n",
-      "text/plain": [
-       "<Figure size 432x288 with 1 Axes>"
-      ]
-     },
-     "metadata": {
-      "needs_background": "light"
-     },
-     "output_type": "display_data"
     }
    ],
    "source": [
-    "bern_drawset.theta.plot.density()"
+    "bern_drawset.shape, bern_drawset.columns"
    ]
   },
   {
@@ -377,16 +441,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 12,
+   "execution_count": 34,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "(4000, 1)"
+       "(5000, 1)"
       ]
      },
-     "execution_count": 12,
+     "execution_count": 34,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -398,7 +462,7 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 13,
+   "execution_count": 35,
    "metadata": {},
    "outputs": [
     {
@@ -428,15 +492,15 @@
        "  <tbody>\n",
        "    <tr>\n",
        "      <th>0</th>\n",
-       "      <td>0.266562</td>\n",
+       "      <td>0.177184</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>1</th>\n",
-       "      <td>0.217474</td>\n",
+       "      <td>0.171893</td>\n",
        "    </tr>\n",
        "    <tr>\n",
        "      <th>2</th>\n",
-       "      <td>0.211124</td>\n",
+       "      <td>0.551869</td>\n",
        "    </tr>\n",
        "  </tbody>\n",
        "</table>\n",
@@ -444,12 +508,12 @@
       ],
       "text/plain": [
        "      theta\n",
-       "0  0.266562\n",
-       "1  0.217474\n",
-       "2  0.211124"
+       "0  0.177184\n",
+       "1  0.171893\n",
+       "2  0.551869"
       ]
      },
-     "execution_count": 13,
+     "execution_count": 35,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -458,6 +522,38 @@
     "thetas[0:3]"
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": 36,
+   "metadata": {},
+   "outputs": [
+    {
+     "data": {
+      "text/plain": [
+       "<matplotlib.axes._subplots.AxesSubplot at 0x1182f1da0>"
+      ]
+     },
+     "execution_count": 36,
+     "metadata": {},
+     "output_type": "execute_result"
+    },
+    {
+     "data": {
+      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYUAAAD8CAYAAACYebj1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4wLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvqOYd8AAAIABJREFUeJzt3Xt83HWd7/HXJ5N7c7+0TZM0aemFFkpbmkK5LchF8QZ6RAUVxUURdddV93LYy3F3dfec1T27e0TcVUQW9KHrBVFR0BUQgVIohFIKvaclbZK2SZrm2tyT7/ljJnHaJs00nd/85vJ+Ph7z6G9mvvnNu8kkn/n+vt/v72fOOURERADS/A4gIiLxQ0VBREQmqSiIiMgkFQUREZmkoiAiIpNUFEREZJKKgoiITFJREBGRSSoKIiIyKd3vAGeqrKzM1dbW+h1DRCShvPzyy0edc+UztUu4olBbW0t9fb3fMUREEoqZHYiknQ4fiYjIJBUFERGZpKIgIiKTVBRERGSSioKIiExSURARkUkqCiIiMinh1ilI6jnSPchjrx0G4NoV81hYmutzIpHkpaIgcW1Tw1Hu+O7L9A2NAvC/H9vJJ65czOeuXUZ6QB1dkWhTUZC49cbR49zx3ZdZUJTNNz60jsz0NL76xF6+/tQ+tjV387Vb1lKUm+l3TJGkoo9aEpecc3zh569jwAMfvYjF5XlUFefyz+9dzZffs4oX9ndwwz3PsetIj99RRZKKioLEpY0NR3l271E+d90yFhTlnPDc+9cv5Ad3XMLgyBjv/vom7nt2P6Nj4z4lFUkuKgoSl7759H7m5mfxwQ0Lp3x+XU0xv/jjy7nknFL+4dGdvP3ujTzXcDTGKUWSj4qCxJ29rb1sbDjKbZfVkpUemLbdvIJsvv2ROr556zr6R0b54H2b+cR36+keGIlhWpHkoqIgceehLc2kpxnvq6uesa2Z8Zbz5vP4567kz9+ynN/uauP933ye3kEVBpHZUFGQuDI6Ns5Pt7Rw1fK5lOVlRfx12RkBPv2mJXzrw3Xsbevjrp+85mFKkeSloiBxZdO+Dtp6h7hpXeWsvv6q5XP5/HXLePS1wzy1qy3K6USSn4qCxJX/3n6E3MwAVy2fO+t9fPyKxSwum8M/PraT8XEXxXQiyU9FQeLG+Ljj8R2tXLW8nOyM6QeYZ5KZnsafXLuUhrY+nlRvQeSMqChI3Hi1uYu23iGuWznvrPf19lUVVBbl8O2N+6OQTCR1qChI3Hh8RyuBNOPq5WdfFNIDadxyUTUv7D9G07H+KKQTSQ0qChI3ntjZysWLSijMzYjK/t61NjhY/dNXWqKyP5FUoKIgceFI9yB7Wvu4cll51PZZVZzLJYtLeXhLM85pwFkkEp4VBTPLNrMXzexVM9tuZn8/RZssM/uhmTWY2WYzq/Uqj8S3iVNUXL60LKr7vWHNAho7+tl1pDeq+xVJVl72FIaAq51zq4E1wPVmtuGkNrcDnc65JcC/AV/2MI/EsY0NRymdk8mK+QVR3e81K4JTW5/Y0RrV/YokK8+KggvqC93NCN1O7sPfCDwY2n4IuMbMzKtMEp+cc2xsOMqlS8pIS4vuj39ufjZrqot4fKeKgkgkPB1TMLOAmW0F2oDHnXObT2pSCTQBOOdGgW6g1MtMEn92t/bS3jvEFUuie+hownUr57GtuZsj3YOe7F8kmXhaFJxzY865NUAVcJGZnT+b/ZjZHWZWb2b17e3t0Q0pvnthXwcAly7x5vPAxLqHJ9RbEJlRTGYfOee6gKeA6096qgWoBjCzdKAQ6Jji6+91ztU55+rKy6M3O0Xiw0sHOllQmE1Vca4n+186N4/Kohye3asPFCIz8XL2UbmZFYW2c4DrgF0nNXsE+Eho+ybgt05zB1OKc46XGztZV1vi2WuYGZcvKWPTvg7GdC4kkdPysqdQATxlZtuAlwiOKfzSzL5oZjeE2nwbKDWzBuDzwF0e5pE41NI1wJGeQdbXFnv6OpctLaN3cJTXWro9fR2RRJfu1Y6dc9uAtVM8/oWw7UHgvV5lkPhX39gJBC+v6aVLzwmOVzzXcJQ11UWevpZIItOKZvFV/YFj5GWlc26U1yecrCwvixUVBbqOs8gMVBTEV/WNnaxdWEQgyusTpnL5klLqGzsZGB7z/LVEEpWKgvhmYHiMPa29MTucc9mSMobHxqk/cCwmryeSiFQUxDc7j/Qw7uD8ysKYvN762hICacbm/SoKItNRURDfbA/NBFoVo6IwJyud8xYU8GKjioLIdFQUxDevtXRTMieTisLsmL3m+toStjZ1MTSqcQWRqagoiG9eb+nhvAUFxPIciOtrSxgeHee1Zq1XEJmKioL4Ymg0OMgcq0NHEyYWyekQksjUVBTEF7uP9DI67mI2yDyhNC+Lc8rn8NIbKgoiU1FREF9sP9QDwPkLYlsUAC5aVEL9gU6dB0lkCioK4os9rb3kZgaoKs6J+Wuvry2hd3CU3bpEp8gpVBTEF3tb+1gyNy/qV1qLxPrQGVlf0riCyClUFMQXe9t6WTo335fXrirOoaIwW4PNIlNQUZCY6x4YobVniKXz8nx5fTNjfW0J9Y3H0OU7RE6koiAx19AWPJa/dK4/RQGCU1Nbe4Zo7hzwLYNIPFJRkJjb29oHwLJ5/hw+AqjTuILIlFQUJOb2tvWRkxGgsij2M48mLJuXT352Oi+FLvIjIkEqChJze1p7fZt5NCGQZqyrKeZlnUZb5AQqChJzDW3B6ah+q6spZk9rH139w35HEYkbKgoSUwPDYxzuHuSc8jl+R5kcV3j5gA4hiUxQUZCYOnisH4CaUv+LwuqqIjICpnEFkTCeFQUzqzazp8xsh5ltN7M/maLNVWbWbWZbQ7cveJVH4kNjx3EAauOgKORkBji/slDjCiJh0j3c9yjwp865LWaWD7xsZo8753ac1O5Z59w7PMwhcaTxaLAoLCzN9TlJUF1NMQ9uOsDgyBjZGQG/44j4zrOegnPusHNuS2i7F9gJVHr1epIYGjv6KZmTSWFOht9RgOC4wvDYOK+36KI7IhCjMQUzqwXWApunePoSM3vVzH5lZufFIo/450DHcWripJcAwZ4CoHEFkRDPi4KZ5QE/AT7rnOs56ektQI1zbjXwNeBn0+zjDjOrN7P69vZ2bwOLpw509MfFeMKE0rwsFpfP0biCSIinRcHMMggWhO855x4++XnnXI9zri+0/RiQYWZlU7S71zlX55yrKy8v9zKyeGhwZIxD3QNx1VMAWF8TvOjOuC66I+Lp7CMDvg3sdM796zRt5ofaYWYXhfJ0eJVJ/NXc2Y9z8THzKNy62mK6+kfY197ndxQR33k5++gy4FbgNTPbGnrsr4CFAM65bwA3AZ80s1FgALjZ6VzGSavx6MQahTjrKUyeHK+TpT6epE8kHnhWFJxzG4HTntzGOXcPcI9XGSS+xNMahXC1pbmU5WVSf+AYH7h4od9xRHylFc0SM40dxynMyaB4TqbfUU5gZtTVlFCvGUgiKgoSOwc6+uPu0NGEutpiDh7rp7Vn0O8oIr5SUZCYae4coLokXotCcFxBvQVJdSoKEhPj446WzgGqfLywzumct6CA7Iw06rVeQVKcioLExNG+IYbHxqkqjs+ikBFIY211sXoKkvJUFCQmmjoHAKiM06IAsL62mO2HuukbGvU7iohvVBQkJlq6gkWhqjg+xxQA1tWWMO5g68Euv6OI+EZFQWKiuTO4cK0yTscUAC5cWESaoXEFSWkqChITLZ0DFOdmMCfLy0X0Zyc/O4Nz5xdoXEFSmoqCxERz50BcjydMWF9bzJaDnYyOjfsdRcQXKgoSEy1dA1QVxe94woR1tSX0D4+x83Cv31FEfKGiIJ5zztHc2Z8QPYWJi+5sOahDSJKaVBTEc8eODzM4Er9rFMItKMqhojCblw+oKEhqUlEQzzVPrFGI45lH4S5cWKyiIClLRUE8lwhrFMJdWFNMS9cAR7p1cjxJPSoK4rnJNQoJcPgIYJ3GFSSFqSiI51o6B8jPTqcwJ8PvKBFZWVFAVnqaDiFJSlJREM81dw4kzHgCQGZ6GquritRTkJSkoiCea+kaSJjxhAkX1hTzeks3gyNjfkcRiSkVBfFUcI3CQEJMRw23rqaYkTHH6y3dfkcRiSkVBfFUz8AofUOjCVcULlxYBKBxBUk5nhUFM6s2s6fMbIeZbTezP5mijZnZ3WbWYGbbzOxCr/KIP5oS4OyoUynNy6KmNJdXdBptSTFe9hRGgT91zq0ENgCfNrOVJ7V5K7A0dLsD+A8P84gPEm2NQrgLqorY1qyiIKnFs6LgnDvsnNsS2u4FdgKVJzW7EfiOC3oBKDKzCq8ySew1J8AV16azuqqQQ92DtPcO+R1FJGZiMqZgZrXAWmDzSU9VAk1h95s5tXBIAmvpHCA3M0BxbmKsUQh3QVVwXOG1FvUWJHV4XhTMLA/4CfBZ51zPLPdxh5nVm1l9e3t7dAOKp5o7+6ksysHM/I5yxs5bUECawatNmoEkqcPTomBmGQQLwveccw9P0aQFqA67XxV67ATOuXudc3XOubry8nJvwoongmsUEu/QEcCcrHSWzM3jNU1LlRTi5ewjA74N7HTO/es0zR4BPhyahbQB6HbOHfYqk8RecI1C4g0yT1hVGRxsds75HUUkJiIqCmb2sJm93czOpIhcBtwKXG1mW0O3t5nZnWZ2Z6jNY8B+oAH4FvCpMwkv8a1vaJTugZGEHGSesLq6kKN9wxzSGVMlRUR6FfV/Bz4K3G1mPwb+0zm3+3Rf4JzbCJz2QLILfvz6dIQZJME0J+gahXCTg83NXQn9/xCJVESf/J1zTzjnPghcCDQCT5jZJjP7aGjcQOQULQk8HXXCufPzSU8zXm3WuIKkhogPB5lZKXAb8DHgFeCrBIvE454kk4T3+4VriVsUsjMCnFuRr0VskjIiOnxkZj8FlgPfBd4ZNhj8QzOr9yqcJLaWzgEy09Mom5Pld5SzsqqyiF9uO4RzLiGn1oqciUh7Ct9yzq10zv2fiYJgZlkAzrk6z9JJQmvuCl5HIS0tsf+Qrq4qpHdwlMaOfr+jiHgu0qLwD1M89nw0g0jySbSL60xnYrBZh5AkFZy2KJjZfDNbB+SY2VozuzB0uwpI3MnnEhMtSVIUls3LIzsjja1NKgqS/GYaU3gLwcHlKiB8AVov8FceZZIkMDgyxtG+oYSeeTQhPZDGqspCFQVJCactCs65B4EHzew9zrmfxCiTJIFDoZlHydBTAFhTXcSDzx9geHSczHRdm0qS10yHjz4U2qw1s8+ffItBPklQyTAdNdya6mKGR8fZdWRW53QUSRgzfeSZE/o3D8if4iYypWRYuBZudXUhAK/qEJIkuZkOH30z9O/fxyaOJIvmzgECacb8gmy/o0RFZVEOZXlZvNLUxa2X+J1GxDuRnhDvK2ZWYGYZZvakmbWHHVoSOUVL1wDzC7JJDyTH8XczY021Bpsl+UX6G/vm0AVy3kHw3EdLgD/3KpQkvmSZjhpuTXUR+9uP0z0w4ncUEc9EWhQmDjO9Hfixc05nB5PTaukaSJrxhAlrqosBLWKT5BZpUfilme0C1gFPmlk5oBPMy5RGx8Y50jOYNDOPJqyqKsQMXjmooiDJK9JTZ98FXArUOedGgOPAjV4Gk8R1pGeQsXGXdIePCnMyWD4vnxffOOZ3FBHPRHqRHYBzCa5XCP+a70Q5jySB5iSbjhpuw+JSfvDSQS1ik6QV6eyj7wL/F7gcWB+66eyoMqXJNQpJ1lOAYFEYHBnXuIIkrUh7CnXASqerl0sEJlYzL0jConDxohIAXtjfQV1tic9pRKIv0v7v68B8L4NI8mjpHKAsL4vsjIDfUaKueE4m587PZ2PDUb+jiHgi0p5CGbDDzF4EhiYedM7d4EkqSWgtXQNJN/Mo3FXL53Lfs/vpHhihMEeXKJfkEmlR+Lsz3bGZ3U9wsVubc+78KZ6/Cvg58EbooYedc18809eR+NPSNcDKBQV+x/DMtSvm8o2n9/HMnnbeuXqB33FEoirSKalPE1zJnBHafgnYMsOXPQBcP0ObZ51za0I3FYQkMD7uaOkcoCoJxxMmrF1YTHFuBk/ubPU7ikjURTr76OPAQ8A3Qw9VAj873dc4554BNKE7xRztG2J4bDwpp6NOCKQZbzp3Lk/tbmdodMzvOCJRFelA86eBy4AeAOfcXmBuFF7/EjN71cx+ZWbnTdfIzO4ws3ozq29vb4/Cy4pXmpPs4jrTeefqBXQPjPDULr0fJblEWhSGnHPDE3dCC9jOdnrqFqDGObca+Bqn6Xk45+51ztU55+rKy8vP8mXFS8l2HYXpXLGkjPL8LB7e0ux3FJGoirQoPG1mfwXkmNl1wI+BX5zNCzvnepxzfaHtx4AMMys7m32K/1pSpKeQHkjjXWsW8NTuNo72Dc38BSIJItKicBfQDrwGfAJ4DPibs3lhM5tvZhbaviiUpeNs9in+a+kcoDAng/zs5J+q+f711YyMOb6/+aDfUUSiJqIpqc65cTP7GfAz51xEB1HN7L+Aq4AyM2sG/hbICO3vG8BNwCfNbBQYAG7WiunE19zZn/S9hAlL5uZz5bJyvvP8AT5x5WKy0pNvsZ6kntMWhdAn+b8F/ohQr8LMxoCvzTSF1Dl3ywzP3wPcc0ZpJe61dA1QUzpn5oZJ4vbLF/Hh+1/kF68e5qZ1VX7HETlrMx0++hzBWUfrnXMlzrkS4GLgMjP7nOfpJKE455Lyimunc8XSMpbNy+PbG99AHV1JBjMVhVuBW5xzE6uOcc7tBz4EfNjLYJJ4ugdGOD48ltSnuDiZmfGHly1i5+EeNus6C5IEZioKGc65U878FRpXSP6RRDkjTceCM4+qinN9ThJb71pbSXFuBvdvfGPmxiJxbqaiMDzL5yQFNXX2A1Bdkjo9BYDsjAAfuHghj+9spelYv99xRM7KTEVhtZn1THHrBVbFIqAkjok/iKnWUwC4dUMtATMe2NTodxSRs3LaouCcCzjnCqa45TvndPhITtDU2U9BdnpKnk56fmE2bzl/Pj/Z0szw6LjfcURmTReZlahp7hyguiT1egkTblpXRVf/CE/tbvM7isisqShI1DQd66c6BQ8dTbhiSRlleTofkiQ2FQWJCudcqKeQWoPM4dIDady4ZgG/3dVGV7/mYUhiUlGQqGjvHWJodDylDx8BvHttJSNjjl+/fsTvKCKzoqIgUTE5HTWFDx8BnLeggOqSHH6zQ1dlk8SkoiBRMbFwLZUPH0FwhfObV85nY8NR+oZG/Y4jcsZUFCQqJtYoVBaldk8B4M0r5zE8Os4ze3RVNkk8KgoSFc2dA5TlZZGTqdNHr6sppmROJr/ZrnEFSTwqChIVTZ39KX/oaEJ6II2rz53Lb3e1MTKmhWySWFQUJCqaOlN7jcLJrl0xl57BUeobO/2OInJGVBTkrI2OjXOoa1A9hTBXLC0nM5DGkzs1C0kSi4qCnLXD3YOMjTv1FMLMyUpnwzmlPLGzVRffkYSioiBnbWKNQiqeHfV0rlsxl8aOfva1H/c7ikjEVBTkrDVrjcKUrl4xD0CHkCShqCjIWWvu7CfNYEEKXZs5EpVFOayoKODJnTprqiQOz4qCmd1vZm1m9vo0z5uZ3W1mDWa2zcwu9CqLeKupc4CKwhwyAvqMcbJrV8yl/sAxOo/rBHmSGLz8LX4AuP40z78VWBq63QH8h4dZxENNx/qpKlYvYSrXrJjHuIPf7VFvQRKDZ0XBOfcMcOw0TW4EvuOCXgCKzKzCqzzineDCNQ0yT+WCykLK87N4QoeQJEH42d+vBJrC7jeHHjuFmd1hZvVmVt/ervPJxJPBkTFae4Y0HXUaaWnG1cvn8szudl2mUxJCQhwEds7d65yrc87VlZeX+x1HwrR0BWce6fDR9K5ZMZfeoVFeajxdx1kkPvhZFFqA6rD7VaHHJIFMnB1Vh4+md/nSMjLT03hCU1MlAfhZFB4BPhyahbQB6HbOHfYxj8xCU6fWKMwkNzOdy84p5cmdbVrdLHHPyymp/wU8Dyw3s2Yzu93M7jSzO0NNHgP2Aw3At4BPeZVFvHOw4ziZ6WnMy8/2O0pcu2bFPA4e66ehrc/vKCKnle7Vjp1zt8zwvAM+7dXrS2wc6OinpiSXtDTzO0pcu2bFXP7mZ/Dr14+wdF6+33FEppUQA80Svw509FNTqvGEmVQU5nDxohJ++kqLDiFJXFNRkFlzznHg2HFqSuf4HSUhvOfCKvYfPc4rTV1+RxGZloqCzFpb7xCDI+PqKUToravmk52RxsNbmv2OIjItFQWZtQMdwemo6ilEJj87g7ecN59Hth5iYHjM7zgiU1JRkFlr7AheJ6BGaxQi9oGLFtIzOMrPt2pJjsQnFQWZtYMd/QTSjEqtZo7YRYtKWFFRwAObGjXgLHFJRUFmrbHjOJVFOmX2mTAzbru0hl1Hetn8hk57IfFHv80yawePaTrqbNy4ppKi3AweeK7R7ygip1BRkFlrPHpcRWEWsjMCfOCihfxmxxEOhgbrReKFioLMSlf/MD2Do9Rq5tGsfOTSWgJpxv3PveF3FJETqCjIrExMR12omUezMq8gmxtWV/Kj+ia6+nWpTokfKgoyKxPTUWvL1FOYrY9dsYj+4TG+t/mg31FEJqkoyKwcVE/hrK2oKOCKpWU8uKlRV2WTuKGiILOy/+hxKgqzyc4I+B0loX38isW09Q7xyKuH/I4iAqgoyCzta+9jydw8v2MkvCuWlnHu/Hzue3a/FrNJXFBRkDPmnGNfWx/nlKsonC0z4/bLF7HrSC/P7j3qdxwRFQU5c609QxwfHuOccg0yR8MNaxZQnp/Ft57d73cUERUFOXP72oOXlFRPITqy0gPcdmktz+49yq4jPX7HkRSnoiBnbLIoaEwhaj548UJyMgLc96wWs4m/VBTkjDW09ZGflc7c/Cy/oySNotxM3ldXxc+3ttDaM+h3HElhnhYFM7vezHabWYOZ3TXF87eZWbuZbQ3dPuZlHomOfe19LJ6bh5n5HSWp/OHlixgddzy4qdHvKJLCPCsKZhYAvg68FVgJ3GJmK6do+kPn3JrQ7T6v8kj07Gs7rkFmD9SUzuEtK+fzvc0H6R8e9TuOpCgvewoXAQ3Ouf3OuWHgB8CNHr6exEDv4AhHegY1yOyR269YRPfACD99RVdmE394WRQqgaaw+82hx072HjPbZmYPmVm1h3kkCva09gKwfF6+z0mSU11NMasqC/nP53RlNvGH3wPNvwBqnXMXAI8DD07VyMzuMLN6M6tvb2+PaUA50c7DwaJwboWKghfMjI9eVktDW58Ws4kvvCwKLUD4J/+q0GOTnHMdzrmh0N37gHVT7cg5d69zrs45V1deXu5JWInMriM95GenU1mk6zJ75e0XVFCen6VrLYgvvCwKLwFLzWyRmWUCNwOPhDcws4qwuzcAOz3MI1Gw63AvK+YXaOaRh7LSA9y6oYbf7W6noa3P7ziSYjwrCs65UeCPgP8m+Mf+R8657Wb2RTO7IdTsM2a23cxeBT4D3OZVHjl7zjl2HenVoaMY+MDFC8kMpGl6qsRcupc7d849Bjx20mNfCNv+S+Avvcwg0dPcOUDf0Cjnzi/wO0rSK8vL4sY1C3jo5Wb+9M3LKMrN9DuSpAi/B5olgew4HDwvj3oKsXH7FYsYGBnj/o0aW5DYUVGQiL3a1EV6mrGyQj2FWDh3fgFvWzWf+59rpPO4ruMssaGiIBF7tbmLcyvydbW1GPrstcs4PjzKvTqttsSIioJEZHzcsa2pmzXVRX5HSSnL5uXzzgsW8MBzjRzp1onyxHsqChKR/Uf76B0aZXWVikKs/dmblzPuHF96dIffUSQFqChIRLY2dQOop+CDhaW5fOqqJTy67TDP7tWKfvGWioJE5OUDx8jPTmexToTni09cuZjFZXP4i4e20dWvQWfxjoqCRGTTvg42LC4lkKaVzH7Izgjw1ZvXcrRviL94aJtOlieeUVGQGTV39nOgo59Lzyn1O0pKW1VVyP+8/lx+s6OVf3tir99xJEl5uqJZksPz+zoAuPScMp+TyO2XL2JPay93P7mXqqIc3rdeZ5uX6FJRkBlt2tdB6ZxMls3TeILfzIx/fPcqDncPctfD28hMT+Nda6e6TInI7OjwkZzWyNg4T+1u4w+WlevMqHEiI5DGN29dx8WLSvncj7by4/qmmb9IJEIqCnJaL75xjK7+Ed5y3ny/o0iY3Mx07r9tPZcvKePPH9rGN57ep8FniQoVBTmtX79+hJyMAFcu08WN4k1OZoBvfbiOd1xQwT/9ahf/6+evMzo27ncsSXAaU5BpjY6N8+vtR7hqeTk5mTrfUTzKzghw981rqSzO4ZtP7+dQ1yD/7+Y1FGRn+B1NEpR6CjKt3+1up713iBvXLPA7ipxGWprxl29dwZfedT5P72nnhq9tZPuhbr9jSYJSUZBpff/Fg5TnZ3HNinl+R5EI3Lqhhh/esYHBkXH+x79v4j+fe4OxcY0zyJlRUZAp7T7Sy293tXHLRQvJCOhtkijqakt49DOXc+k5pfz9L3bwrq8/x3MNRzUILRHTmIJM6e7f7mVOZoCPXlrrdxQ5Q6V5Wdx/23p+ue0w//DoDj5432ZWVhTwjtUVXLyohBUVBeRm6ldfpqZ3hpzihf0dPLrtMJ+5egnFc3Rt4ERkZrxz9QKuWzmPn2xp5kf1zXzl17snn59fkM3C0lxqSnKpKc1lYekcakpyWVw+h3wNUqc0S7RuZV1dnauvr/c7RtLqGRzhnV/byNi44/HPXalZR0mktWeQbc3d7DzcQ2PHcQ529HPgWD/tvUOTbdLTjA2LS7lu5TyuWzmPBUU5PiaWaDKzl51zdTO287IomNn1wFeBAHCfc+6fTno+C/gOsA7oAN7vnGs83T5VFLwzODLGx79Tz/P7OvjhJzawrqbE70gSA/3Doxw8Fjzp4SsHu3h8xxH2tR8H4LwFBVy7IlggVlYUkKaz5CYs34uCmQWAPcB1QDPwEnCLc25HWJtPARc45+40s5uBdzvn3n+6/aooeKPpWD+f+cErvHKwi6/cdAHvq9OJ1lLZvvY+Ht/RyhOZzbEnAAAH7klEQVQ7Wnn5YCfOQX5WOhdUF7Kqsoja0lwWluRSVZxLSV4mczIDOg1KnIu0KHg5pnAR0OCc2x8K9APgRiD8moI3An8X2n4IuMfMzCXaMa0ENDo2zuHuQXYc7uHxHa08svUQWelpfP0DF/L2Cyr8jic+O6c8j3OuzOPOK8/haN8QT+9uZ8vBTrY2dXHfs/sZPWmqa0bAKM7NpGROJkW5GRTmZJCXlUFeVoC87PTgdnY6+Vnp5GWlhx5LZ05WOtkZaWSnB8jOCJCVnqbeiM+8LAqVQPiZupqBi6dr45wbNbNuoBQ4Gu0wT+9p50u/DNajiZoz+bYOe39PbIbXpd8/NnE/7Dl34r/hTnmdmfYxRbuJR92ZZpziufAX6B8Zm5zDnpeVznvrqvj0m5boGLKcoiwvi/esq+I966qA33+gaDrWT0vXAJ39wxw7PkJX/zDHjg/T2T/MgY5+egdH6RsK3s5kvURmII2sjDSyMwJkZ6QRCPVAJnoikyXDTvjn1OeT0PvXV/OxKxZ7+hoJMfvIzO4A7gBYuHDhrPaRl5XO8nn5YTs94Z8Tur6/f+yU5lO/8Sb3FbaPU/Yf3txOfezkd/ps9jHF/+XUfQU3cjMDVBXnsLg8jzXVRWSmay2CRCY9kEZ1SS7VJbkRtXfOMTgyPlkg+gZH6R0aoW9wlOPDowyOjDM0Msbg6DiDI2MMjgT/HRoNbo87d8oHp+k+2DkiLz6JqCwvy/PX8LIotADhB6arQo9N1abZzNKBQoIDzidwzt0L3AvBMYXZhFlXU8y6muLZfKmInAUzIyczQE5mgPJ87/+oydnx8uPhS8BSM1tkZpnAzcAjJ7V5BPhIaPsm4LcaTxAR8Y9nPYXQGMEfAf9NcErq/c657Wb2RaDeOfcI8G3gu2bWABwjWDhERMQnno4pOOceAx476bEvhG0PAu/1MoOIiEROo4siIjJJRUFERCapKIiIyCQVBRERmaSiICIikxLu1Nlm1g4cOMMvK8ODU2dEQbzmAmWbrXjNFq+5QNlmYza5apxz5TM1SriiMBtmVh/J2QFjLV5zgbLNVrxmi9dcoGyz4WUuHT4SEZFJKgoiIjIpVYrCvX4HmEa85gJlm614zRavuUDZZsOzXCkxpiAiIpFJlZ6CiIhEICmLgpmVmNnjZrY39O+0F1IwswIzazaze+Ihl5mtMbPnzWy7mW0zs9NeszoKma43s91m1mBmd03xfJaZ/TD0/GYzq/Uyzxlm+7yZ7Qh9n540s5p4yBXW7j1m5swsZrNXIslmZu8Lfd+2m9n34yWbmS00s6fM7JXQz/RtMcp1v5m1mdnr0zxvZnZ3KPc2M7swFrkizPbBUKbXzGyTma0+6xd1ziXdDfgKcFdo+y7gy6dp+1Xg+8A98ZALWAYsDW0vAA4DRR7lCQD7gMVAJvAqsPKkNp8CvhHavhn4YYx+hpFkexOQG9r+ZCyyRZIr1C4feAZ4AaiLo+/ZUuAVoDh0f24cZbsX+GRoeyXQGKNsfwBcCLw+zfNvA35F8OKGG4DNscgVYbZLw36Wb41GtqTsKQA3Ag+Gth8E3jVVIzNbB8wDfhMvuZxze5xze0Pbh4A2YMYFJ7N0EdDgnNvvnBsGfhDKOF3mh4BrbKrrffqQzTn3lHOuP3T3BYJX9/M9V8iXgC8DgzHIdCbZPg583TnXCeCca4ujbA4oCG0XAodiEcw59wzB67lM50bgOy7oBaDIzCriIZtzbtPEz5Io/Q4ka1GY55w7HNo+QvAP/wnMLA34F+DP4ilXODO7iOCnqn0e5akEmsLuN4cem7KNc24U6AZKPcpzptnC3U7w05zXZswVOrxQ7Zx7NAZ5wkXyPVsGLDOz58zsBTO7Po6y/R3wITNrJngdlj+OTbQZnel70S9R+R3w9CI7XjKzJ4D5Uzz11+F3nHPOzKaaYvUp4DHnXHM0P/hGIdfEfiqA7wIfcc6NRy1gEjKzDwF1wJVxkCUN+FfgNp+jTCed4CGkqwh+qnzGzFY557p8TRV0C/CAc+5fzOwSgldlPF/v/5mZ2ZsIFoXLz3ZfCVsUnHPXTvecmbWaWYVz7nDoj+tUXeRLgCvM7FNAHpBpZn3OuWkHDmOUCzMrAB4F/jrUXfVKC1Addr8q9NhUbZrNLJ1gt77Dw0xnkg0zu5Zgwb3SOTcUB7nygfOB34U+bMwHHjGzG5xz9T5ng+Cn3M3OuRHgDTPbQ7BIvBQH2W4Hrgdwzj1vZtkEz/ETq0Nc04novegXM7sAuA94q3Pu7H83YzVgEssb8M+cOKD7lRna30ZsBppnzEXwcNGTwGdjkCcd2A8s4veDf+ed1ObTnDjQ/KMY/QwjybaW4KG1pTF8b82Y66T2vyN2A82RfM+uBx4MbZcRPCxSGifZfgXcFtpeQXBMwWL0vatl+sHct3PiQPOLsXq/RZBtIdAAXBq114vlfy6G38TS0B/WvcATQEno8Trgvinax6oozJgL+BAwAmwNu63xMNPbgD2hP65/HXrsi8ANoe1s4MehN96LwOIY/hxnyvYE0Br2fXokHnKd1DZmRSHC75kRPLy1A3gNuDmOsq0EngsVjK3Am2OU678IzvIbIdiTuh24E7gz7Hv29VDu12L885wp231AZ9jvQP3ZvqZWNIuIyKRknX0kIiKzoKIgIiKTVBRERGSSioKIiExSURARkUkqCiIiMklFQUREJqkoiIjIpP8P2Zb+NYi/cZoAAAAASUVORK5CYII=\n",
+      "text/plain": [
+       "<Figure size 432x288 with 1 Axes>"
+      ]
+     },
+     "metadata": {
+      "needs_background": "light"
+     },
+     "output_type": "display_data"
+    }
+   ],
+   "source": [
+    "bern_drawset.theta.plot.density()"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -473,16 +569,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 14,
+   "execution_count": 37,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "(1000, 4, 8)"
+       "(1000, 5, 8)"
       ]
      },
-     "execution_count": 14,
+     "execution_count": 37,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -492,18 +588,30 @@
     "bern_fit.sample.shape"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Python's index slicing operations can be used to access the information by chain.\n",
+    "For example, to select all draws and all output columns from the first chain,\n",
+    "we specify the chain index (2nd index dimension).  As arrays indexing starts at 0,\n",
+    "the index '0' corresponds to the first chain in the ``RunSet``.\n",
+    "\n",
+    "The following expression selects the first 3 draws from chain 1 for the parameter `theta`:"
+   ]
+  },
   {
    "cell_type": "code",
-   "execution_count": 15,
+   "execution_count": 38,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "('theta', array([0.266562, 0.217474, 0.211124]))"
+       "('theta', array([0.177184, 0.171893, 0.551869]))"
       ]
      },
-     "execution_count": 15,
+     "execution_count": 38,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -525,16 +633,16 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 16,
+   "execution_count": 39,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "array([0.750967, 0.962835, 1.09088 , 0.947622])"
+       "array([1.39734 , 1.27207 , 0.956467, 0.804156, 0.867051])"
       ]
      },
-     "execution_count": 16,
+     "execution_count": 39,
      "metadata": {},
      "output_type": "execute_result"
     }
@@ -556,19 +664,20 @@
   },
   {
    "cell_type": "code",
-   "execution_count": 17,
+   "execution_count": 40,
    "metadata": {},
    "outputs": [
     {
      "data": {
       "text/plain": [
-       "('diag_e', array([[0.66346 ],\n",
-       "        [0.507504],\n",
-       "        [0.493591],\n",
-       "        [0.427833]]))"
+       "('diag_e', array([[0.444845],\n",
+       "        [0.441921],\n",
+       "        [0.609734],\n",
+       "        [0.610979],\n",
+       "        [0.566994]]))"
       ]
      },
-     "execution_count": 17,
+     "execution_count": 40,
      "metadata": {},
      "output_type": "execute_result"
     }
diff --git a/docs/_static/css/project-template.css b/docs/_static/css/project-template.css
index f6caff2..479a349 100644
--- a/docs/_static/css/project-template.css
+++ b/docs/_static/css/project-template.css
@@ -1,5 +1,9 @@
 @import url("theme.css");
 
+.highlight * {
+    background: #eff2f5;
+}
+
 .highlight a {
     text-decoration: underline;
 }
diff --git a/docs/api.rst b/docs/api.rst
index 12d3dc6..9159616 100644
--- a/docs/api.rst
+++ b/docs/api.rst
@@ -1,6 +1,6 @@
-###
-API
-###
+#############
+API Reference
+#############
 
 ********
 Commands
@@ -11,42 +11,46 @@ Model fit and analysis
 
 compile_model
 -------------
-.. autofunction:: cmdstanpy.cmds.compile_model
+.. autofunction:: cmdstanpy.compile_model
 
 sample
 ------
-.. autofunction:: cmdstanpy.cmds.sample
+.. autofunction:: cmdstanpy.sample
 
 get_drawset
 -----------
-.. autofunction:: cmdstanpy.cmds.get_drawset
+.. autofunction:: cmdstanpy.get_drawset
 
 save_csvfiles
 -------------
-.. autofunction:: cmdstanpy.cmds.save_csvfiles
+.. autofunction:: cmdstanpy.save_csvfiles
 
 summary
 -------
-.. autofunction:: cmdstanpy.cmds.summary
+.. autofunction:: cmdstanpy.summary
 
 diagnose
 --------
-.. autofunction:: cmdstanpy.cmds.diagnose
+.. autofunction:: cmdstanpy.diagnose
 
 Utility functions
 =================
 
 set_cmdstan_path
 ----------------
-.. autofunction:: cmdstanpy.utils.set_cmdstan_path
+.. autofunction:: cmdstanpy.set_cmdstan_path
 
 cmdstan_path
 ------------
-.. autofunction:: cmdstanpy.utils.cmdstan_path
+.. autofunction:: cmdstanpy.cmdstan_path
 
 jsondump
 --------
-.. autofunction:: cmdstanpy.utils.jsondump
+.. autofunction:: cmdstanpy.jsondump
+
+rdump
+--------
+.. autofunction:: cmdstanpy.rdump
 
 
 *******
@@ -57,11 +61,11 @@ Classes
 Model
 =====
 
-.. autoclass:: cmdstanpy.lib.Model
+.. autoclass:: cmdstanpy.Model
    :members:
 
 RunSet
 ======
 
-.. autoclass:: cmdstanpy.lib.RunSet
+.. autoclass:: cmdstanpy.RunSet
    :members:
diff --git a/docs/conf.py b/docs/conf.py
index 43f69f3..9b20cf8 100644
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -82,7 +82,7 @@ copyright = '2019, Stan Development Team'
 # built documents.
 #
 # The short X.Y version.
-from cmdstanpy._version import __version__
+from cmdstanpy import __version__
 version = __version__
 # The full version, including alpha/beta/rc tags.
 release = __version__
diff --git a/docs/getting_started.rst b/docs/getting_started.rst
new file mode 100644
index 0000000..794b056
--- /dev/null
+++ b/docs/getting_started.rst
@@ -0,0 +1,175 @@
+Getting Started
+===============
+
+
+Installation
+____________
+
+CmdStanPy is a pure-Python package which can be installed from PyPI
+
+.. code-block:: bash
+
+    pip install --upgrade cmdstanpy
+
+or from GitHub
+
+.. code-block:: bash
+
+    pip install -e git+https://github.com/stan-dev/cmdstanpy
+
+CmdStanPy requires a local install of CmdStan.
+If you don't have CmdStan installed, you can run the CmdStanPy script ``install_cmdstan``
+which downloads CmdStan from GitHub and builds the CmdStan utilities.
+By default this script installs the latest version of CmdStan into a directory named
+``.cmdstanpy`` in your ``$HOME`` directory:
+
+.. code-block:: bash
+
+    install_cmdstan
+    ls -F ~/.cmdstanpy
+
+The named arguments: `-d <directory>` and  `-v <version>`
+can be used to override these defaults:
+
+.. code-block:: bash
+
+    install_cmdstan -d my_local_cmdstan -v 2.19.1
+    ls -F my_local_cmdstan
+
+If you already have CmdStan installed in a directory
+then you can set the environment variable ``CMDSTAN`` to this
+location and it will be picked up by CmdStanPy:
+
+.. code-block:: bash
+
+    export CMDSTAN='/path/to/cmdstan-2.19.1'
+
+
+The CmdStanPy commands ``cmdstan_path`` and ``set_cmdstan_path``
+get and set this environment variable:
+
+.. code-block:: python
+
+    from cmdstanpy import cmdstan_path, set_cmdstan_path
+
+    oldpath = cmdstan_path()
+    set_cmdstan_path(os.join('path','to','cmdstan'))
+    newpath = cmdstan_path()
+
+
+
+CmdStanPy's "Hello World"
+_________________________
+
+To exercise the essential functions of CmdStanPy, we will
+compile and run the example Stan model ``bernoulli.stan`` which is
+distributed with CmdStan.
+
+
+Specify a Stan model
+--------------------
+
+The ``Model`` class specifies the Stan program and its corresponding compiled executable.
+The function ``compile_model`` is used to compile or or recompile a Stan program.
+It takes the path to the Stan program file and returns a ``Model`` object:
+
+.. code-block:: python
+
+    import os
+    import os.path
+    from cmdstanpy import Model, compile_model
+
+    bernoulli_stan = os.path.join(cmdstan_path(), 'examples', 'bernoulli', 'bernoulli.stan')
+    bernoulli_model = compile_model(bernoulli_stan)
+
+If you already have a compiled executable, you can construct a ``Model`` object directly:
+
+.. code-block:: python
+
+    bernoulli_model = Model(
+            stan_file=os.path.join(cmdstan_path(), 'examples', 'bernoulli', 'bernoulli.stan')
+            stan_exe=os.path.join(cmdstan_path(), 'examples', 'bernoulli', 'bernoulli')
+            )
+
+
+Run the HMC-NUTS sampler
+------------------------
+
+The ``sample`` function invokes the Stan HMC-NUTS sampler on the ``Model`` object and some data
+and returns a ``RunSet`` object:
+
+.. code-block:: python
+
+    bern_data = { "N" : 10, "y" : [0,1,0,0,0,0,0,0,0,1] }
+    bern_fit = sample(bernoulli_model, data=bern_data)
+    
+By default, the ``sample`` command runs 4 sampler chains.
+The ``RunSet`` object records the results of each sampler chain.
+If no output file path is specified, the sampler outputs
+are written to a temporary directory which is deleted
+when the current Python session is terminated.
+
+
+Summarize or save the results
+-----------------------------
+
+The ``get_drawset`` function is used to get the draws from
+all chains as a ``pandas.DataFrame``, one draw per row, one column per
+model parameter, transformed parameter, generated quantity variable.
+The ``params`` argument is used to restrict the DataFrame
+columns to just the specified parameter names.
+
+.. code-block:: python
+
+    get_drawset(bern_fit, params=['theta'])
+
+Underlyingly, this information is stored in the ``sample`` property
+of a ``RunSet`` object as a 3-D ``numpy.ndarray`` (i.e., a multi-dimensional array)
+with dimensions: (draws, chains, columns).
+Python's index slicing operations can be used to access the information by chain.
+For example, to select all draws and all output columns from the first chain,
+we specify the chain index (2nd index dimension).  As arrays indexing starts at 0,
+the index '0' corresponds to the first chain in the ``RunSet``:
+
+.. code-block:: python
+
+    chain_1 = bern_fit.sample[:,0,:]
+
+
+CmdStan is distributed with a posterior analysis utility ``stansummary``
+that reads the outputs of all chains and computes summary statistics
+on the model fit for all parameters. CmdStanPy's ``summary`` function
+runs the CmdStan ``stansummary`` utility and returns the output as a pandas.DataFrame:
+
+.. code-block:: python
+
+    summary(bern_fit)
+
+CmdStan is distributed with a second posterior analysis utility ``diagnose``
+that reads the outputs of all chains and checks for the following
+potential problems:
+
++ Transitions that hit the maximum treedepth
++ Divergent transitions
++ Low E-BFMI values (sampler transitions HMC potential energy)
++ Low effective sample sizes
++ High R-hat values
+
+The ``diagnose`` function runs the CmdStan ``diagnose`` utility
+and prints the output to the console.
+
+.. code-block:: python
+
+    diagnose(bern_fit)
+
+By default, CmdStanPy will save all CmdStan outputs in a temporary
+directory which is deleted when the Python session exits.
+In particular, unless the ``csv_output_file`` argument to the ``sample``
+function is overtly specified, all the csv output files will be written into
+this temporary directory and then when the session exits.
+The ``save_csvfiles`` function moves the CmdStan csv output files
+to the specified location, renaming them using a specified basename.
+
+.. code-block:: python
+
+    save_csvfiles(bern_fit, dir='some/path', basename='descriptive-name')
diff --git a/docs/guide.rst b/docs/guide.rst
deleted file mode 100644
index af448e1..0000000
--- a/docs/guide.rst
+++ /dev/null
@@ -1,135 +0,0 @@
-Guide
-=====
-
-This tutorial will walk you through the different steps of using CmdStanPy. 
-
-Install
-_______
-
-CmdStanPy can be installed from GitHub
-
-.. code-block:: bash
-
-	pip install -e git+https://github.com/stan-dev/cmdstanpy
-
-CmdStanPy requires a local install of CmdStan.
-If you don't have CmdStan installed, you can run the script ``make_cmdstan.sh`` which
-will download CmdStan from GitHub and build the CmdStan utilities.
-By default this script installs the latest version of CmdStan into a directory named
-`.cmdstanpy` in the user's `$HOME` directory:
-
-
-.. code-block:: bash
-
-	./make_cmdstan.sh
-    ls -F ~/.cmdstanpy
-
-This script takes two optional named arguments: `-d <directory> -v <version>`
-
-.. code-block:: bash
-
-	./make_cmdstan.sh -d cmdstan -v 2.18.1
-    ls -F ~/cmdstan
-
-If you already have CmdStan installed, then set the ``CMDSTAN`` environment variable accordingly,
-either from within your Python session or setting the environment variable directly using ``bash``:
-
-.. code-block:: bash
-
-	export CMDSTAN='/path/to/cmdstan'
-
-.. code-block:: python
-
-	import os
-	os.environ['CMDSTAN'] = '/path/to/cmdstan'
-
-
-Basics
-______
-
-The ``compile_model`` function takes as its argument the name of the Stan program and returns a ``Model`` object:
-
-.. code-block:: python
-
-	import os
-	import os.path
-	from cmdstanpy.lib import Model
-	from cmdstanpy.cmds import compile_model
-	bernoulli_stan = os.path.join('cmdstanpy', 'test', 'data', 'bernoulli.stan')
-	bernoulli_model = compile_model(bernoulli_stan)
-	print(bernoulli_model)
-	bernoulli_model.name
-
-The ``Model`` class specifies the Stan program and its corresponding compiled executable.
-If you already have a compiled executable, you can construct the Model object directly:
-
-.. code-block:: python
-
-	bernoulli_model = Model(
-			stan_file=os.path.join('cmdstanpy', 'test', 'data', 'bernoulli.stan'),
-			stan_exe=os.path.join('cmdstanpy', 'test', 'data', 'bernoulli')
-		    )
-	print(bernoulli_model)
-	bernoulli_model.name
-
-The ``sample`` function invokes the Stan HMC-NUTS sampler on the ``Model`` object and some data
-and returns a ``RunSet`` object:
-
-.. code-block:: python
-
-    bern_data = { "N" : 10, "y" : [0,1,0,0,0,0,0,0,0,1] }
-    bern_fit = sample(bernoulli_model, chains=4, cores=2, data=bern_data)
-
-The ``sample`` property of the ``RunSet`` object is a 3-D ``numpy.ndarray``
-which contains all draws across all chains, stored column major format so that values
-for each parameter are stored contiguously in memory.
-The dimensions of the ndarray are arranged (draws, chains, columns).
-
-The ``get_drawset`` function flattens this 3-D ndarray to a pandas.DataFrame,
-one draw per row.  The `params` argument is used to restrict the DataFrame
-view to the specified parameter names, else all output columns are returned.
-
-.. code-block:: python
-
-    bern_fit.sample.shape
-    get_drawset(bern_fit, params=['theta'])
-
-
-CmdStan is distributed with a posterior analysis utility `stansummary`
-that reads the outputs of all chains and computes summary statistics
-on the model fit for all parameters. CmdStanPy's ``summary`` function
-runs this utility and returns the output as a pandas.DataFrame:
-
-.. code-block:: python
-
-    summary(bern_fit)
-
-CmdStan is distributed with a second posterior analysis utility `diagnose`
-that reads the outputs of all chains and checks for the following
-potential problems:
-
-+ Transitions that hit the maximum treedepth
-+ Divergent transitions
-+ Low E-BFMI values (sampler transitions HMC potential energy)
-+ Low effective sample sizes
-+ High R-hat values
-
-The ``diagnose`` function prints the output of the CmdStan ``bin/diagnose``:
-
-.. code-block:: python
-
-    diagnose(bern_fit)
-
-By default, CmdStanPy will save all CmdStan outputs in a temporary
-directory which is deleted when the Python session exits.
-In particular, if the ``sample`` command is invoked without
-specifying the `csv_output_file` path, then the csv output files
-will be written into this temporary directory and therefore will
-be deleted once the session exits.
-The ``save_csvfiles`` function moves the CmdStan csv output files
-to the specified location, renaming them using a specified basename.
-
-.. code-block:: python
-
-    save_csvfiles(bern_fit, dir='some/path', basename='descriptive-name')
-
diff --git a/docs/index.rst b/docs/index.rst
index af4ecfb..fd32f76 100644
--- a/docs/index.rst
+++ b/docs/index.rst
@@ -14,15 +14,12 @@
 
 
 .. toctree::
-   :maxdepth: 3
+   :maxdepth: 4
    :caption: Contents:
 
-   guide
+   getting_started
    api
 
 
-Indices and tables
-==================
+:ref:`genindex`
 
-* :ref:`genindex`
-* :ref:`search`
diff --git a/make_cmdstan.sh b/make_cmdstan.sh
deleted file mode 100755
index 2b84ec9..0000000
--- a/make_cmdstan.sh
+++ /dev/null
@@ -1,76 +0,0 @@
-#!/usr/bin/env bash
-
-# install a cmdstan release into a specified directory
-#  - build binaries, compile example model to build model header
-#  - symlink downloaded version as "cmdstan"
-
-while getopts ":d:v:" opt; do
-  case $opt in
-    d) RELDIR="$OPTARG"
-    ;;
-    v) VER="$OPTARG"
-    ;;
-    \?) echo "Invalid option -$OPTARG" >&2
-    ;;
-  esac
-done
-
-if [ -z ${RELDIR} ]; then
-    RELDIR="$HOME/.cmdstanpy"
-fi
-
-if [[ ! -e ${RELDIR} ]]; then
-   mkdir ${RELDIR}
-fi
-if [[ ! -d ${RELDIR} ]]; then
-    echo "cannot install cmdstan, ${RELDIR} is not a directory"
-    exit 1
-fi
-echo "cmdstan dir: ${RELDIR}"
-
-if [ -z ${VER} ]; then
-    TAG=`curl -s https://api.github.com/repos/stan-dev/cmdstan/releases/latest | grep "tag_name"`
-    echo $TAG > tmp-tag
-    VER=`perl -p -e 's/"tag_name": "v//g; s/",//g' tmp-tag`
-    rm tmp-tag
-fi
-CS=cmdstan-${VER}
-echo "cmdstan version: ${VER}"
-
-pushd ${RELDIR} > /dev/null
-if [[ -d $cs && -f ${CS}/bin/stanc && -f ${CS}/examples/bernoulli/bernoulli ]]; then
-    echo "cmdstan already installed"
-    exit 0
-fi
-echo "installing ${CS}"
-
-echo "downloading ${CS}.tar.gz from github"
-curl -s -OL https://github.com/stan-dev/cmdstan/releases/download/v${VER}/${CS}.tar.gz
-CURL_RC=$?
-if [[ ${CURL_RC} -ne 0 ]]; then
-    echo "github download failed, curl exited with: ${CURL_RC}"
-    exit ${CURL_RC}
-fi
-echo "download complete"
-
-echo "unpacking archive"
-tar xzf ${CS}.tar.gz
-TAR_RC=$?
-if [[ ${TAR_RC} -ne 0 ]]; then
-    echo "corrupt download file ${CS}.tar.gz, tar exited with: ${TAR_RC}"
-    exit ${TAR_RC}
-fi
-
-if [[ -h cmdstan ]]; then
-    unlink cmdstan
-fi
-ln -s ${CS} cmdstan
-pushd cmdstan > /dev/null
-echo "building cmdstan binaries"
-make -j2 build examples/bernoulli/bernoulli
-echo "installed ${CS}; symlink: `ls -l ${RELDIR}/cmdstan`"
-
-# cleanup
-pushd -0 > /dev/null
-dirs -c > /dev/null
-echo "all installed versions in ${RELDIR} dir: `ls -Fd ${RELDIR}/* | grep -v @ | grep -v gz`"
diff --git a/setup.py b/setup.py
index 7ab2669..33d99b3 100644
--- a/setup.py
+++ b/setup.py
@@ -1,7 +1,7 @@
 #!/usr/bin/env python
 
 import setuptools
-
+from cmdstanpy import __version__
 
 def readme_contents() -> str:
     with open('README.md', 'r') as fd:
@@ -13,7 +13,7 @@ _classifiers = """
 Programming Language :: Python :: 3
 License :: OSI Approved :: Apache Software License
 Operating System :: OS Independent
-Development Status :: Alpha
+Development Status :: 4 - Beta
 Intended Audience :: Science/Research
 Natural Language :: English
 Programming Language :: Python
@@ -37,13 +37,14 @@ EXTRAS_REQUIRE = {
 
 setuptools.setup(
     name='cmdstanpy',
-    version='0.10',
+    version=__version__,
     description='Python interface to CmdStan',
     long_description=readme_contents(),
     long_description_content_type="text/markdown",
     author='Stan Dev Team',
     url='https://github.com/stan-dev/cmdstanpy',
     packages=['cmdstanpy'],
+    scripts=['bin/install_cmdstan'],
     install_requires=INSTALL_REQUIRES,
     extras_require=EXTRAS_REQUIRE,
     classifiers=_classifiers.strip().split('\n'),


Test patch:
diff --git a/test/test_utils.py b/test/test_utils.py
index c142dcc..b5870dd 100644
--- a/test/test_utils.py
+++ b/test/test_utils.py
@@ -5,7 +5,7 @@ import json
 
 from cmdstanpy import TMPDIR
 from cmdstanpy.lib import StanData
-from cmdstanpy.utils import cmdstan_path, set_cmdstan_path, validate_cmdstan_path, check_csv
+from cmdstanpy.utils import cmdstan_path, set_cmdstan_path, validate_cmdstan_path, get_latest_cmdstan, check_csv
 
 
 datafiles_path = os.path.join('test', 'data')
@@ -18,19 +18,21 @@ y <- c(0, 1, 0, 0, 0, 0, 0, 0, 0, 1)
 class CmdStanPathTest(unittest.TestCase):
     def test_default_path(self):
         abs_rel_path = os.path.expanduser(os.path.join('~', '.cmdstanpy', 'cmdstan'))
-        self.assertEqual(abs_rel_path, cmdstan_path())
+        self.assertTrue(cmdstan_path().startswith(abs_rel_path))
 
     def test_set_path(self):
-        abs_rel_path = os.path.expanduser(os.path.join('~', '.cmdstanpy', 'cmdstan'))
-        self.assertEqual(abs_rel_path, cmdstan_path())
+        install_dir = os.path.expanduser(os.path.join('~', '.cmdstanpy'))
         install_version = os.path.expanduser(
-            os.path.join('~', '.cmdstanpy', 'cmdstan-2.19.1'))
+            os.path.join(install_dir, get_latest_cmdstan(install_dir)))
         set_cmdstan_path(install_version)
         self.assertEqual(install_version, cmdstan_path())
 
     def test_validate_path(self):
-        abs_rel_path = os.path.expanduser(os.path.join('~', '.cmdstanpy', 'cmdstan'))
-        validate_cmdstan_path(abs_rel_path)
+        install_dir = os.path.expanduser(os.path.join('~', '.cmdstanpy'))
+        install_version = os.path.expanduser(
+            os.path.join(install_dir, get_latest_cmdstan(install_dir)))
+        set_cmdstan_path(install_version)
+        validate_cmdstan_path(install_version)
         path_foo = os.path.abspath(os.path.join('releases', 'foo'))
         with self.assertRaisesRegex(ValueError, 'no such CmdStan directory'):
             validate_cmdstan_path(path_foo)