补充-其他一些模块介绍

# Seaborn¶

Matplotlib has proven to be an incredibly useful and popular visualization tool, but even avid users will admit it often leaves much to be desired. There are several valid complaints about Matplotlib that often come up:

- Prior to version 2.0, Matplotlib's defaults are not exactly the best choices. It was based off of MATLAB circa 1999, and this often shows.
- Matplotlib's API is relatively low level. Doing sophisticated statistical visualization is possible, but often requires a lot of boilerplate code.
- Matplotlib predated Pandas by more than a decade, and thus is not designed for use with Pandas `DataFrame`s. In order to visualize data from a Pandas `DataFrame`, you must extract each `Series` and often concatenate them together into the right format. It would be nicer to have a plotting library that can intelligently use the `DataFrame` labels in a plot.

An answer to these problems is Seaborn. Seaborn provides an API on top of Matplotlib that offers sane choices for plot style and color defaults, defines simple high-level functions for common statistical plot types, and integrates with the functionality provided by Pandas `DataFrame`s.

## Seaborn Versus Matplotlib¶

Here is an example of a simple random-walk plot in Matplotlib, using its classic plot formatting and colors. We start with the typical imports:

In [1]:

    import matplotlib.pyplot as plt
    plt.style.use('classic')
    %matplotlib inline
    import numpy as np
    import pandas as pd

    # Create some data
    rng = np.random.RandomState(0)
    x = np.linspace(0, 10, 500)
    y = np.cumsum(rng.randn(500, 6), 0)

    # Plot the data with Matplotlib defaults
    plt.plot(x, y)
    plt.legend('ABCDEF', ncol=2, loc='upper left');

In [2]:

    import seaborn as sns
    sns.set()

In [3]:

    # same plotting code as above!
    plt.plot(x, y)
    plt.legend('ABCDEF', ncol=2, loc='upper left');

## Histograms, KDE, and densities¶

Often in statistical data visualization, all you want is to plot histograms and joint distributions of variables. We have seen that this is relatively straightforward in Matplotlib:

In [4]:

    data = np.random.multivariate_normal([0, 0], [[5, 2], [2, 2]], size=2000)
    data = pd.DataFrame(data, columns=['x', 'y'])

    for col in 'xy':
        plt.hist(data[col], normed=True, alpha=0.5)

In [5]:

    for col in 'xy':
        sns.kdeplot(data[col], shade=True)

    /Users/lee/anaconda/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

In [6]:

    sns.distplot(data['x'])
    sns.distplot(data['y']);

    /Users/lee/anaconda/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

In [7]:

    sns.kdeplot(data);

In [8]:

    with sns.axes_style('white'):
        sns.jointplot("x", "y", data, kind='kde');

    /Users/lee/anaconda/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

In [9]:

    with sns.axes_style('white'):
        sns.jointplot("x", "y", data, kind='hex')

## Pair plots¶

When you generalize joint plots to datasets of larger dimensions, you end up with pair plots. This is very useful for exploring correlations between multidimensional data, when you'd like to plot all pairs of values against each other.

We'll demo this with the well-known Iris dataset, which lists measurements of petals and sepals of three iris species:

In [10]:

    iris = sns.load_dataset("iris")
    iris.head()

Out[10]:

| sepal_length | sepal_width | petal_length | petal_width | species
--- | --- | --- | --- | --- | ---
0 | 5.1 | 3.5 | 1.4 | 0.2 | setosa
1 | 4.9 | 3.0 | 1.4 | 0.2 | setosa
2 | 4.7 | 3.2 | 1.3 | 0.2 | setosa
3 | 4.6 | 3.1 | 1.5 | 0.2 | setosa
4 | 5.0 | 3.6 | 1.4 | 0.2 | setosa

In [11]:

    sns.pairplot(iris, hue='species', size=2.5);

## Faceted histograms¶

In [12]:

    tips = sns.load_dataset('tips')
    tips.head()

Out[12]:

| total_bill | tip | sex | smoker | day | time | size
--- | --- | --- | --- | --- | --- | --- | ---
0 | 16.99 | 1.01 | Female | No | Sun | Dinner | 2
1 | 10.34 | 1.66 | Male | No | Sun | Dinner | 3
2 | 21.01 | 3.50 | Male | No | Sun | Dinner | 3
3 | 23.68 | 3.31 | Male | No | Sun | Dinner | 2
4 | 24.59 | 3.61 | Female | No | Sun | Dinner | 4

In [13]:

    tips['tip_pct'] = 100 * tips['tip'] / tips['total_bill']

    grid = sns.FacetGrid(tips, row="sex", col="time", margin_titles=True)
    grid.map(plt.hist, "tip_pct", bins=np.linspace(0, 40, 15));

## Factor plots¶

Factor plots can be useful for this kind of visualization as well. This allows you to view the distribution of a parameter within bins defined by any other parameter:

In [14]:

    with sns.axes_style(style='ticks'):
        g = sns.factorplot("day", "total_bill", "sex", data=tips, kind="box")
        g.set_axis_labels("Day", "Total Bill");

## Joint distributions¶

Similar to the pairplot we saw earlier, we can use `sns.jointplot` to show the joint distribution between different datasets, along with the associated marginal distributions:

In [15]:

    with sns.axes_style('white'):
        sns.jointplot("total_bill", "tip", data=tips, kind='hex')

The joint plot can even do some automatic kernel density estimation and regression:

In [16]:

    sns.jointplot("total_bill", "tip", data=tips, kind='reg');

    /Users/lee/anaconda/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

## Bar plots¶

Time series can be plotted using `sns.factorplot`.

In [17]:

    planets = sns.load_dataset('planets')
    planets.head()

Out[17]:

| method | number | orbital_period | mass | distance | year
--- | --- | --- | --- | --- | --- | ---
0 | Radial Velocity | 1 | 269.300 | 7.10 | 77.40 | 2006
1 | Radial Velocity | 1 | 874.774 | 2.21 | 56.95 | 2008
2 | Radial Velocity | 1 | 763.000 | 2.60 | 19.84 | 2011
3 | Radial Velocity | 1 | 326.030 | 19.40 | 110.62 | 2007
4 | Radial Velocity | 1 | 516.220 | 10.50 | 119.47 | 2009

In [18]:

    with sns.axes_style('white'):
        g = sns.factorplot("year", data=planets, aspect=2,
                           kind="count", color='steelblue')
        g.set_xticklabels(step=5)

In [19]:

    with sns.axes_style('white'):
        g = sns.factorplot("year", data=planets, aspect=4.0, kind='count',
                           hue='method', order=range(2001, 2015))
        g.set_ylabels('Number of Planets Discovered')

## Example: Exploring Marathon Finishing Times¶

Here we'll look at using Seaborn to help visualize and understand finishing results from a marathon. I've scraped the data from sources on the Web, aggregated it and removed any identifying information, and put it on GitHub where it can be downloaded (if you are interested in using Python for web scraping, I would recommend Web Scraping with Python by Ryan Mitchell). We will start by downloading the data from the Web, and loading it into Pandas:

In [23]:

    # !curl -O https://raw.githubusercontent.com/jakevdp/marathon-data/master/marathon-data.csv

In [25]:

    data = pd.read_csv('https://raw.githubusercontent.com/jakevdp/marathon-data/master/marathon-data.csv')
    data.head()

Out[25]:

| age | gender | split | final
--- | --- | --- | --- | ---
0 | 33 | M | 01:05:38 | 02:08:51
1 | 32 | M | 01:06:26 | 02:09:28
2 | 31 | M | 01:06:49 | 02:10:42
3 | 38 | M | 01:06:16 | 02:13:45
4 | 31 | M | 01:06:32 | 02:13:59

In [26]:

    data.dtypes

Out[26]:

    age        int64
    gender    object
    split     object
    final     object
    dtype: object

In [28]:

    # fix this by providing a converter for the times:
    def convert_time(s):
        h, m, s = map(int, s.split(':'))
        return pd.datetools.timedelta(hours=h, minutes=m, seconds=s)

    data = pd.read_csv('https://raw.githubusercontent.com/jakevdp/marathon-data/master/marathon-data.csv',
                       converters={'split':convert_time, 'final':convert_time})
    data.head()

Out[28]:

| age | gender | split | final
--- | --- | --- | --- | ---
0 | 33 | M | 01:05:38 | 02:08:51
1 | 32 | M | 01:06:26 | 02:09:28
2 | 31 | M | 01:06:49 | 02:10:42
3 | 38 | M | 01:06:16 | 02:13:45
4 | 31 | M | 01:06:32 | 02:13:59

In [29]:

    data.dtypes

Out[29]:

    age                 int64
    gender             object
    split     timedelta64[ns]
    final     timedelta64[ns]
    dtype: object

In [30]:

    data['split_sec'] = data['split'].astype(int) / 1E9
    data['final_sec'] = data['final'].astype(int) / 1E9
    data.head()

Out[30]:

| age | gender | split | final | split_sec | final_sec
--- | --- | --- | --- | --- | --- | ---
0 | 33 | M | 01:05:38 | 02:08:51 | 3938.0 | 7731.0
1 | 32 | M | 01:06:26 | 02:09:28 | 3986.0 | 7768.0
2 | 31 | M | 01:06:49 | 02:10:42 | 4009.0 | 7842.0
3 | 38 | M | 01:06:16 | 02:13:45 | 3976.0 | 8025.0
4 | 31 | M | 01:06:32 | 02:13:59 | 3992.0 | 8039.0

In [31]:

    with sns.axes_style('white'):
        g = sns.jointplot("split_sec", "final_sec", data, kind='hex')
        g.ax_joint.plot(np.linspace(4000, 16000),
                        np.linspace(8000, 32000), ':k')

In [32]:

    data['split_frac'] = 1 - 2 * data['split_sec'] / data['final_sec']
    data.head()

Out[32]:

| age | gender | split | final | split_sec | final_sec | split_frac
--- | --- | --- | --- | --- | --- | --- | ---
0 | 33 | M | 01:05:38 | 02:08:51 | 3938.0 | 7731.0 | -0.018756
1 | 32 | M | 01:06:26 | 02:09:28 | 3986.0 | 7768.0 | -0.026262
2 | 31 | M | 01:06:49 | 02:10:42 | 4009.0 | 7842.0 | -0.022443
3 | 38 | M | 01:06:16 | 02:13:45 | 3976.0 | 8025.0 | 0.009097
4 | 31 | M | 01:06:32 | 02:13:59 | 3992.0 | 8039.0 | 0.006842

Where this split difference is less than zero, the person negative-split the race by that fraction. Let's do a distribution plot of this split fraction:

In [33]:

    sns.distplot(data['split_frac'], kde=False);
    plt.axvline(0, color="k", linestyle="--");

In [34]:

    sum(data.split_frac < 0)

Out[34]:

    251

Out of nearly 40,000 participants, there were only 250 people who negative-split their marathon.

Let's see whether there is any correlation between this split fraction and other variables. We'll do this using a `pairgrid`, which draws plots of all these correlations:

In [35]:

    g = sns.PairGrid(data, vars=['age', 'split_sec', 'final_sec', 'split_frac'],
                     hue='gender', palette='RdBu_r')
    g.map(plt.scatter, alpha=0.8)
    g.add_legend();

The difference between men and women here is interesting. Let's look at the histogram of split fractions for these two groups:

In [36]:

    sns.kdeplot(data.split_frac[data.gender=='M'], label='men', shade=True)
    sns.kdeplot(data.split_frac[data.gender=='W'], label='women', shade=True)
    plt.xlabel('split_frac');

    /Users/lee/anaconda/lib/python3.5/site-packages/statsmodels/nonparametric/kdetools.py:20: VisibleDeprecationWarning: using a non-integer number instead of an integer will result in an error in the future
      y = X[:m/2+1] + np.r_[0,X[m/2+1:],0]*1j

A nice way to compare distributions is to use a violin plot

In [37]:

    sns.violinplot("gender", "split_frac", data=data,
                   palette=["lightblue", "lightpink"]);

In [38]:

    data['age_dec'] = data.age.map(lambda age: 10 * (age // 10))
    data.head()

Out[38]:

| age | gender | split | final | split_sec | final_sec | split_frac | age_dec
--- | --- | --- | --- | --- | --- | --- | --- | ---
0 | 33 | M | 01:05:38 | 02:08:51 | 3938.0 | 7731.0 | -0.018756 | 30
1 | 32 | M | 01:06:26 | 02:09:28 | 3986.0 | 7768.0 | -0.026262 | 30
2 | 31 | M | 01:06:49 | 02:10:42 | 4009.0 | 7842.0 | -0.022443 | 30
3 | 38 | M | 01:06:16 | 02:13:45 | 3976.0 | 8025.0 | 0.009097 | 30
4 | 31 | M | 01:06:32 | 02:13:59 | 3992.0 | 8039.0 | 0.006842 | 30

In [39]:

    men = (data.gender == 'M')
    women = (data.gender == 'W')

    with sns.axes_style(style=None):
        sns.violinplot("age_dec", "split_frac", hue="gender", data=data,
                       split=True, inner="quartile",
                       palette=["lightblue", "lightpink"]);

The 80-year-old women seem to outperform everyone in terms of their split time. This is probably due to the fact that we're estimating the distribution from small numbers, as there are only a handful of runners in that range:

In [40]:

    (data.age > 80).sum()

Out[40]:

    7

Back to the men with negative splits: who are these runners? Does this split fraction correlate with finishing quickly? We can plot this very easily. We'll use `regplot`, which will automatically fit a linear regression to the data:

In [41]:

    g = sns.lmplot('final_sec', 'split_frac', col='gender', data=data,
                   markers=".", scatter_kws=dict(color='c'))
    g.map(plt.axhline, y=0.1, color="k", ls=":");

# Bokeh¶

Bokeh is a Python interactive visualization library that targets modern web browsers for presentation. Bokeh provides elegant, concise construction of novel graphics with high-performance interactivity over very large or streaming datasets in a quick and easy way.

In [1]:

    import bokeh
    bokeh.__version__

Out[1]:

    '0.12.2'

In [2]:

    from bokeh.plotting import figure, output_file, show

    # prepare some data
    x = [1, 2, 3, 4, 5]
    y = [6, 7, 2, 4, 5]

    # output to static HTML file
    output_file("lines.html")

    # create a new plot with a title and axis labels
    p = figure(title="simple line example", x_axis_label='x', y_axis_label='y')

    # add a line renderer with legend and line thickness
    p.line(x, y, legend="Temp.", line_width=2)

    # show the results
    show(p)

In [3]:

    import numpy as np

    from bokeh.plotting import figure, output_file, show

    # prepare some data
    N = 4000
    x = np.random.random(size=N) * 100
    y = np.random.random(size=N) * 100
    radii = np.random.random(size=N) * 1.5
    colors = [
        "#%02x%02x%02x" % (int(r), int(g), 150) for r, g in zip(50+2*x, 30+2*y)
    ]

    # output to static HTML file (with CDN resources)
    output_file("color_scatter.html", title="color_scatter.py example", mode="cdn")

    TOOLS="resize,crosshair,pan,wheel_zoom,box_zoom,reset,box_select,lasso_select"

    # create a new plot with the tools above, and explicit ranges
    p = figure(tools=TOOLS, x_range=(0,100), y_range=(0,100))

    # add a circle renderer with vectorized colors and sizes
    p.circle(x,y, radius=radii, fill_color=colors, fill_alpha=0.6, line_color=None)

    # show the results
    show(p)

In [4]:

    import numpy as np

    from bokeh.layouts import gridplot
    from bokeh.plotting import figure, output_file, show

    # prepare some data
    N = 100
    x = np.linspace(0, 4*np.pi, N)
    y0 = np.sin(x)
    y1 = np.cos(x)
    y2 = np.sin(x) + np.cos(x)

    # output to static HTML file
    output_file("linked_panning.html")

    # create a new plot
    s1 = figure(width=250, plot_height=250, title=None)
    s1.circle(x, y0, size=10, color="navy", alpha=0.5)

    # NEW: create a new plot and share both ranges
    s2 = figure(width=250, height=250, x_range=s1.x_range, y_range=s1.y_range, title=None)
    s2.triangle(x, y1, size=10, color="firebrick", alpha=0.5)

    # NEW: create a new plot and share only one range
    s3 = figure(width=250, height=250, x_range=s1.x_range, title=None)
    s3.square(x, y2, size=10, color="olive", alpha=0.5)

    # NEW: put the subplots in a gridplot
    p = gridplot([[s1, s2, s3]], toolbar_location=None)

    # show the results
    show(p)

In [6]:

    bokeh.sampledata.download()

    Creating /Users/lee/.bokeh directory
    Creating /Users/lee/.bokeh/data directory
    Using data directory: /Users/lee/.bokeh/data
    Downloading: CGM.csv (1589982 bytes)
       1589982 [100.00%]
    Downloading: US_Counties.zip (3182088 bytes)
       3182088 [100.00%]
    Unpacking: US_Counties.csv
    Downloading: us_cities.json (713565 bytes)
        713565 [100.00%]
    Downloading: unemployment09.csv (253301 bytes)
        253301 [100.00%]
    Downloading: AAPL.csv (166698 bytes)
        166698 [100.00%]
    Downloading: FB.csv (9706 bytes)
          9706 [100.00%]
    Downloading: GOOG.csv (113894 bytes)
        113894 [100.00%]
    Downloading: IBM.csv (165625 bytes)
        165625 [100.00%]
    Downloading: MSFT.csv (161614 bytes)
        161614 [100.00%]
    Downloading: WPP2012_SA_DB03_POPULATION_QUINQUENNIAL.zip (5148539 bytes)
       5148539 [100.00%]
    Unpacking: WPP2012_SA_DB03_POPULATION_QUINQUENNIAL.csv
    Downloading: gapminder_fertility.csv (64346 bytes)
         64346 [100.00%]
    Downloading: gapminder_population.csv (94509 bytes)
         94509 [100.00%]
    Downloading: gapminder_life_expectancy.csv (73243 bytes)
         73243 [100.00%]
    Downloading: gapminder_regions.csv (7781 bytes)
          7781 [100.00%]
    Downloading: world_cities.zip (646858 bytes)
        646858 [100.00%]
    Unpacking: world_cities.csv
    Downloading: airports.json (6373 bytes)
          6373 [100.00%]
    Downloading: movies.db.zip (5067833 bytes)
       5067833 [100.00%]
    Unpacking: movies.db

In [7]:

    import numpy as np

    from bokeh.plotting import figure, output_file, show
    from bokeh.sampledata.stocks import AAPL

    # prepare some data
    aapl = np.array(AAPL['adj_close'])
    aapl_dates = np.array(AAPL['date'], dtype=np.datetime64)

    window_size = 30
    window = np.ones(window_size)/float(window_size)
    aapl_avg = np.convolve(aapl, window, 'same')

    # output to static HTML file
    output_file("stocks.html", title="stocks.py example")

    # create a new plot with a a datetime axis type
    p = figure(width=800, height=350, x_axis_type="datetime")

    # add renderers
    p.circle(aapl_dates, aapl, size=4, color='darkgrey', alpha=0.2, legend='close')
    p.line(aapl_dates, aapl_avg, color='navy', legend='avg')

    # NEW: customize by setting attributes
    p.title.text = "AAPL One-Month Average"
    p.legend.location = "top_left"
    p.grid.grid_line_alpha=0
    p.xaxis.axis_label = 'Date'
    p.yaxis.axis_label = 'Price'
    p.ygrid.band_fill_color="olive"
    p.ygrid.band_fill_alpha = 0.1

    # show the results
    show(p)

# Pymc3¶

PyMC3 is a Python package for Bayesian statistical modeling and Probabilistic Machine Learning which focuses on advanced Markov chain Monte Carlo and variational fitting algorithms. Its flexibility and extensibility make it applicable to a large suite of problems.

    pip install git+https://github.com/pymc-devs/pymc3

In [9]:

    import pymc3
    pymc3.__version__

Out[9]:

    '3.1'

## A Motivating Example: Linear Regression¶

To introduce model definition, fitting and posterior analysis, we first consider a simple Bayesian linear regression model with normal priors for the parameters. We are interested in predicting outcomes $Y$ as normally-distributed observations with an expected value $\mu$ that is a linear function of two predictor variables, $X_1$ and $X_2$.

$$Y \sim N(\mu,{\sigma}^2)$$$$\mu = \alpha + {{\beta}_1}X_1 + {{\beta}_2}X_2$$

where $\alpha$ is the intercept, and ${\beta}_i$ is the coefficient for covariate $X_i$, while σσ represents the observation error. Since we are constructing a Bayesian model, the unknown variables in the model must be assigned a prior distribution. We choose zero-mean normal priors with variance of 100 for both regression coefficients, which corresponds to weak information regarding the true parameter values. We choose a half-normal distribution (normal distribution bounded at zero) as the prior for σσ.

$$\alpha \sim N(0,100)$$$${\beta}_i \sim N(0,100) $$$$ \sigma \sim |N(0,1)|$$
In [10]:

    import numpy as np
    import matplotlib.pyplot as plt

    # Initialize random number generator
    np.random.seed(123)

    # True parameter values
    alpha, sigma = 1, 1
    beta = [1, 2.5]

    # Size of dataset
    size = 100

    # Predictor variable
    X1 = np.random.randn(size)
    X2 = np.random.randn(size) * 0.2

    # Simulate outcome variable
    Y = alpha + beta[0]*X1 + beta[1]*X2 + np.random.randn(size)*sigma

In [11]:

    %matplotlib inline

    fig, axes = plt.subplots(1, 2, sharex=True, figsize=(10,4))
    axes[0].scatter(X1, Y)
    axes[1].scatter(X2, Y)
    axes[0].set_ylabel('Y'); axes[0].set_xlabel('X1'); axes[1].set_xlabel('X2');

In [12]:

    import pymc3 as pm

In [13]:

    basic_model = pm.Model()

    with basic_model:

        # Priors for unknown model parameters
        alpha = pm.Normal('alpha', mu=0, sd=10)
        beta = pm.Normal('beta', mu=0, sd=10, shape=2)
        sigma = pm.HalfNormal('sigma', sd=1)

        # Expected value of outcome
        mu = alpha + beta[0]*X1 + beta[1]*X2

        # Likelihood (sampling distribution) of observations
        Y_obs = pm.Normal('Y_obs', mu=mu, sd=sigma, observed=Y)

In [15]:

    help(pm.Normal) #try help(Model), help(Uniform) or help(basic_model)

    Help on class Normal in module pymc3.distributions.continuous:

    class Normal(pymc3.distributions.distribution.Continuous)
     |  Univariate normal log-likelihood.
     |
     |  .. math::
     |
     |     f(x \mid \mu, \tau) =
     |         \sqrt{\frac{\tau}{2\pi}}
     |         \exp\left\{ -\frac{\tau}{2} (x-\mu)^2 \right\}
     |
     |  ========  ==========================================
     |  Support   :math:`x \in \mathbb{R}`
     |  Mean      :math:`\mu`
     |  Variance  :math:`\dfrac{1}{\tau}` or :math:`\sigma^2`
     |  ========  ==========================================
     |
     |  Normal distribution can be parameterized either in terms of precision
     |  or standard deviation. The link between the two parametrizations is
     |  given by
     |
     |  .. math::
     |
     |     \tau = \dfrac{1}{\sigma^2}
     |
     |  Parameters
     |  ----------
     |  mu : float
     |      Mean.
     |  sd : float
     |      Standard deviation (sd > 0).
     |  tau : float
     |      Precision (tau > 0).
     |
     |  Method resolution order:
     |      Normal
     |      pymc3.distributions.distribution.Continuous
     |      pymc3.distributions.distribution.Distribution
     |      builtins.object
     |
     |  Methods defined here:
     |
     |  __init__(self, mu=0, sd=None, tau=None, **kwargs)
     |      Initialize self.  See help(type(self)) for accurate signature.
     |
     |  logp(self, value)
     |
     |  random(self, point=None, size=None, repeat=None)
     |
     |  ----------------------------------------------------------------------
     |  Methods inherited from pymc3.distributions.distribution.Distribution:
     |
     |  __getnewargs__(self)
     |
     |  default(self)
     |
     |  get_test_val(self, val, defaults)
     |
     |  getattr_value(self, val)
     |
     |  ----------------------------------------------------------------------
     |  Class methods inherited from pymc3.distributions.distribution.Distribution:
     |
     |  dist(*args, **kwargs) from builtins.type
     |
     |  ----------------------------------------------------------------------
     |  Static methods inherited from pymc3.distributions.distribution.Distribution:
     |
     |  __new__(cls, name, *args, **kwargs)
     |      Create and return a new object.  See help(type) for accurate signature.
     |
     |  ----------------------------------------------------------------------
     |  Data descriptors inherited from pymc3.distributions.distribution.Distribution:
     |
     |  __dict__
     |      dictionary for instance variables (if defined)
     |
     |  __weakref__
     |      list of weak references to the object (if defined)

In [16]:

    map_estimate = pm.find_MAP(model=basic_model)

    map_estimate

    Optimization terminated successfully.
             Current function value: 149.017982
             Iterations: 16
             Function evaluations: 21
             Gradient evaluations: 21

Out[16]:

    {'alpha': array(0.9065985497559482),
     'beta': array([ 0.94848602,  2.60705514]),
     'sigma_log__': array(-0.032781470174030644)}

In [17]:

    from scipy import optimize

    map_estimate = pm.find_MAP(model=basic_model, fmin=optimize.fmin_powell)

    map_estimate

    Optimization terminated successfully.
             Current function value: 149.019762
             Iterations: 4
             Function evaluations: 176

Out[17]:

    {'alpha': array(0.9090521898977764),
     'beta': array([ 0.95140146,  2.61437458]),
     'sigma_log__': array(-0.030009775203258385)}

In [18]:

    from scipy import optimize

    with basic_model:
        # draw 500 posterior samples
        trace = pm.sample()

    INFO:pymc3:Auto-assigning NUTS sampler...
    INFO:pymc3:Initializing NUTS using ADVI...
    Average Loss = 156.6:   5%|▌         | 10565/200000 [00:01<00:23, 8236.20it/s]
    INFO:pymc3.variational.inference:Convergence archived at 11100
    INFO:pymc3.variational.inference:Interrupted at 11,100 [5%]: Average Loss = 237.04
    100%|██████████| 1000/1000 [00:01<00:00, 855.36it/s]

In [19]:

    _ = pm.traceplot(trace)

In [20]:

    pm.summary(trace)

    alpha:

      Mean             SD               MC Error         95% HPD interval
      -------------------------------------------------------------------

      0.904            0.103            0.004            [0.714, 1.118]

      Posterior quantiles:
      2.5            25             50             75             97.5
      |--------------|==============|==============|--------------|

      0.689          0.839          0.902          0.971          1.114

    beta:

      Mean             SD               MC Error         95% HPD interval
      -------------------------------------------------------------------

      0.949            0.085            0.003            [0.779, 1.117]
      2.613            0.539            0.030            [1.579, 3.643]

      Posterior quantiles:
      2.5            25             50             75             97.5
      |--------------|==============|==============|--------------|

      0.783          0.890          0.948          1.002          1.124
      1.502          2.235          2.648          2.963          3.614

    sigma:

      Mean             SD               MC Error         95% HPD interval
      -------------------------------------------------------------------

      0.992            0.075            0.003            [0.852, 1.133]

      Posterior quantiles:
      2.5            25             50             75             97.5
      |--------------|==============|==============|--------------|

      0.860          0.939          0.987          1.038          1.164

In [ ]: