Merge branch 'master' into version_increment_3.2_alpha

Author: fonnesbeck

.github/ISSUE_TEMPLATE.md

@@ -0,0 +1,19 @@
+If you have questions about a specific use case, please post it to our discourse channel: https://discourse.pymc.io  
+If you are not sure whether this is a bug or not, please also post it on discourse.
+
+## Description of your problem
+
+**Please provide a minimal, self-contained, and reproducible example.**
+
+**Please provide the full traceback.**
+
+**Please provide any additional information below.**
+
+
+## Versions and main components
+
+* PyMC3 Version:
+* Theano Version:
+* Python Version:
+* Operating system:
+* How did you install PyMC3: (conda/pip)

.pylintrc

@@ -47,6 +47,7 @@ enable=import-error,
        used-before-assignment,
        cell-var-from-loop,
        global-variable-undefined,
+       dangerous-default-value,
        # redefined-builtin,
        redefine-in-handler,
        unused-import,

.travis.yml

@@ -30,6 +30,7 @@ env:
   - PYTHON_VERSION=3.6 FLOATX='float64' TESTCMD="--durations=10 --cov-append pymc3/tests/test_variational_inference.py pymc3/tests/test_updates.py"
 script:
   - . ./scripts/test.sh $TESTCMD
+  - . ./scripts/confirm_mpl_optional.sh
 
 after_success:
   - coveralls

README.rst

@@ -10,8 +10,7 @@ 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.
 
-Check out the `getting started
-guide <http://pymc-devs.github.io/pymc3/notebooks/getting_started.html>`__!
+Check out the :ref:`getting started guide<notebooks/getting_started.ipynb>`!
 
 Features
 ========
@@ -38,10 +37,9 @@ Getting started
 If you already know about Bayesian statistics:
 ----------------------------------------------
 
--  `API quickstart guide <http://pymc-devs.github.io/pymc3/notebooks/api_quickstart.html>`__
--  The `PyMC3 tutorial <http://pymc-devs.github.io/pymc3/notebooks/getting_started.html>`__
--  `PyMC3 examples <http://pymc-devs.github.io/pymc3/examples.html>`__
-   and the `API reference <http://pymc-devs.github.io/pymc3/api.html>`__
+-  :ref:`API quickstart guide<notebooks/api_quickstart.ipynb>`
+-  The :ref:`PyMC3 tutorial<notebooks/getting_started.ipynb>`
+-  :ref:`PyMC3 examples<examples>` and the :ref:`API reference<api>`
 
 Learn Bayesian statistics with a book together with PyMC3:
 ----------------------------------------------------------

docs/source/advanced_theano.rst

@@ -99,3 +99,119 @@ Good reasons for defining a custom Op might be the following:
 
 Theano has extensive `documentation, <http://deeplearning.net/software/theano/extending/index.html>`_
 about how to write new Ops.
+
+
+Finding the root of a function
+------------------------------
+
+We'll use finding the root of a function as a simple example.
+Let's say we want to define a model where a parameter is defined
+implicitly as the root of a function, that depends on another
+parameter:
+
+.. math::
+
+   \theta \sim N^+(0, 1)\\
+   \text{$\mu\in \mathbb{R}^+$ such that $R(\mu, \theta)
+         = \mu + \mu e^{\theta \mu} - 1= 0$}\\
+   y \sim N(\mu, 0.1^2)
+
+First, we observe that this problem is well-defined, because
+:math:`R(\cdot, \theta)` is monotone and has the image :math:`(-1, \infty)`
+for :math:`\mu, \theta \in \mathbb{R}^+`. To avoid overflows in
+:math:`\exp(\mu \theta)` for large
+values of :math:`\mu\theta` we instead find the root of
+
+.. math::
+
+    R'(\mu, \theta)
+        = \log(R(\mu, \theta) + 1)
+        = \log(\mu) + \log(1 + e^{\theta\mu}).
+
+Also, we have
+
+.. math::
+
+    \frac{\partial}{\partial\mu}R'(\mu, \theta)
+        = \theta\, \text{logit}^{-1}(\theta\mu) + \mu^{-1}.
+
+We can now use `scipy.optimize.newton` to find the root::
+
+    from scipy import optimize, special
+    import numpy as np
+
+    def func(mu, theta):
+        thetamu = theta * mu
+        value = np.log(mu) + np.logaddexp(0, thetamu)
+        return value
+
+    def jac(mu, theta):
+        thetamu = theta * mu
+        jac = theta * special.expit(thetamu) + 1 / mu
+        return jac
+
+    def mu_from_theta(theta):
+        return optimize.newton(func, 1, fprime=jac, args=(0.4,))
+
+We could wrap `mu_from_theta` with `tt.as_op` and use gradient-free
+methods like Metropolis, but to get NUTS and ADVI working, we also
+need to define the derivative of `mu_from_theta`. We can find this
+derivative using the implicit function theorem, or equivalently we
+take the derivative with respect of :math:`\theta` for both sides of
+:math:`R(\mu(\theta), \theta) = 0` and solve for :math:`\frac{d\mu}{d\theta}`.
+This isn't hard to do by hand, but for the fun of it, let's do it using
+sympy::
+
+    import sympy
+
+    mu = sympy.Function('mu')
+    theta = sympy.Symbol('theta')
+    R = mu(theta) + mu(theta) * sympy.exp(theta * mu(theta)) - 1
+    solution = sympy.solve(R.diff(theta), mu(theta).diff(theta))[0]
+
+We get
+
+.. math::
+
+    \frac{d}{d\theta}\mu(\theta)
+        = - \frac{\mu(\theta)^2}{1 + \theta\mu(\theta) + e^{-\theta\mu(\theta)}}
+
+Now, we use this to define a theano op, that also computes the gradient::
+
+    import theano
+    import theano.tensor as tt
+    import theano.tests.unittest_tools
+
+    class MuFromTheta(tt.Op):
+        itypes = [tt.dscalar]
+        otypes = [tt.dscalar]
+
+        def perform(self, node, inputs, outputs):
+            theta, = inputs
+            mu = mu_from_theta(theta)
+            outputs[0][0] = np.array(mu)
+
+        def grad(self, inputs, g):
+            theta, = inputs
+            mu = self(theta)
+            thetamu = theta * mu
+            return [- g[0] * mu ** 2 / (1 + thetamu + tt.exp(-thetamu))]
+
+If you value your sanity, always check that the gradient is ok::
+
+    theano.tests.unittest_tools.verify_grad(MuFromTheta(), [np.array(0.2)])
+    theano.tests.unittest_tools.verify_grad(MuFromTheta(), [np.array(1e-5)])
+    theano.tests.unittest_tools.verify_grad(MuFromTheta(), [np.array(1e5)])
+
+We can now define our model using this new op::
+
+    import pymc3 as pm
+
+    tt_mu_from_theta = MuFromTheta()
+
+    with pm.Model() as model:
+        theta = pm.HalfNormal('theta', sd=1)
+        mu = pm.Deterministic('mu', tt_mu_from_theta(theta))
+        pm.Normal('y', mu=mu, sd=0.1, observed=[0.2, 0.21, 0.3])
+
+        trace = pm.sample()

docs/source/api/gp.rst

@@ -1,28 +1,35 @@
-******************
-Gaussian Processes
-******************
-
-
-GP
---
+GP Implementations
+------------------
 
 .. currentmodule:: pymc3.gp.gp
 .. automodule:: pymc3.gp.gp
    :members:
 
+Mean Functions
+--------------
+
+.. currentmodule:: pymc3.gp.cov
+.. autosummary::
+
+   Zero
+   Constant
+   Linear
 
-Covariance Functions / Kernels
-------------------------------
+Covariance Functions
+--------------------
 
 .. currentmodule:: pymc3.gp.cov
 .. autosummary::
 
+   Constant
+   WhiteNoise
    ExpQuad
    RatQuad
    Matern32
    Matern52
    Exponential
    Cosine
+   Periodic
    Linear
    Polynomial
    WarpedInput

docs/source/api/plots.rst

@@ -5,4 +5,5 @@ Plots
 .. currentmodule:: pymc3.plots
 
 .. automodule:: pymc3.plots
-   :members: traceplot, plot_posterior, forestplot, compare_plot, autocorrplot
+   :members: traceplot, plot_posterior, forestplot, compareplot, autocorrplot,
+             energyplot, kdeplot

docs/source/examples.rst

@@ -13,6 +13,7 @@ Howto
    notebooks/posterior_predictive.ipynb
    notebooks/model_comparison.ipynb
    notebooks/model_averaging.ipynb
+   notebooks/Bayes_factor.ipynb
    notebooks/howto_debugging.ipynb
    notebooks/PyMC3_tips_and_heuristic.ipynb
    notebooks/LKJ.ipynb
@@ -48,8 +49,8 @@ Gaussian Processes
 ==================
 
 .. toctree::
+   notebooks/GP-MeansAndCovs.ipynb
    notebooks/GP-Marginal.ipynb
-   notebooks/GP-Covariances.ipynb
    notebooks/GP-Latent.ipynb
    notebooks/GP-SparseApprox.ipynb
    notebooks/GP-TProcess.ipynb

docs/source/getting_started.rst

@@ -8,5 +8,6 @@ Getting started
    notebooks/getting_started.ipynb
    notebooks/api_quickstart.ipynb
    notebooks/variational_api_quickstart.ipynb
+   gp.rst
    theano.rst
    advanced_theano.rst

docs/source/gp.rst

@@ -1,6 +1,6 @@
-===========================
-Gaussian Processes in PyMC3
-===========================
+******************
+Gaussian Processes
+******************
 
 GP Basics
 =========
@@ -111,8 +111,8 @@ with non-standard covariance and mean functons.  For more information check out
 the tutorial on covariance functions.
 
 
-:code:`gp.*` implementations
-============================
+GP Implementations
+==================
 
 PyMC3 includes several GP implementations, including marginal and latent
 variable models and also some fast approximations.  Their usage all follows a