Skip to content

Commit

Permalink
removing extra assert in tests
Browse files Browse the repository at this point in the history
  • Loading branch information
dfm committed Oct 19, 2017
1 parent 750a82d commit ff4e5ec
Show file tree
Hide file tree
Showing 5 changed files with 91 additions and 319 deletions.
243 changes: 55 additions & 188 deletions docs/_static/notebooks/quickstart.ipynb

Large diffs are not rendered by default.

1 change: 1 addition & 0 deletions docs/index.rst
Original file line number Diff line number Diff line change
Expand Up @@ -73,6 +73,7 @@ If you run into any issues, please don't hesitate to `open an issue on GitHub
:caption: User Guide

user/install
user/faq

.. toctree::
:maxdepth: 1
Expand Down
113 changes: 32 additions & 81 deletions docs/tutorials/quickstart.rst
Original file line number Diff line number Diff line change
Expand Up @@ -22,7 +22,7 @@ This notebook was made with the following version of emcee:
.. parsed-literal::
'0.3.0.dev0'
'3.0.0.dev0'
Expand Down Expand Up @@ -60,9 +60,9 @@ requires the logarithm of :math:`p`. We’ll call it ``log_prob``:

.. code:: python
def log_prob(x, mu, icov):
def log_prob(x, mu, cov):
diff = x - mu
return -0.5*np.dot(diff,np.dot(icov,diff))
return -0.5*np.dot(diff, np.linalg.solve(cov,diff))
It is important that the first argument of the probability function is
the position of a single "walker" (a *N* dimensional ``numpy`` array).
Expand All @@ -85,24 +85,17 @@ dimensions:
cov += cov.T - np.diag(cov.diagonal())
cov = np.dot(cov,cov)
and where ``cov`` is :math:`\Sigma`. Before going on, let's compute the
inverse of ``cov`` because that's what we need in our probability
function:
and where ``cov`` is :math:`\Sigma`.

.. code:: python
icov = np.linalg.inv(cov)
It's probably overkill this time but how about we use 250 walkers?
Before we go on, we need to guess a starting point for each of the 250
walkers. This position will be a 50-dimensional vector so the initial
guess should be a 250-by-50 array—or a list of 250 arrays that each have
50 elements. It's not a very good guess but we'll just guess a random
number between 0 and 1 for each component:
How about we use 32 walkers? Before we go on, we need to guess a
starting point for each of the 32 walkers. This position will be a
5-dimensional vector so the initial guess should be a 32-by-5 array.
It's not a very good guess but we'll just guess a random number between
0 and 1 for each component:

.. code:: python
nwalkers = 250
nwalkers = 32
p0 = np.random.rand(nwalkers, ndim)
Now that we've gotten past all the bookkeeping stuff, we can move on to
Expand All @@ -111,15 +104,15 @@ the fun stuff. The main interface provided by ``emcee`` is the

.. code:: python
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=[means, icov])
sampler = emcee.EnsembleSampler(nwalkers, ndim, log_prob, args=[means, cov])
Remember how our function ``log_prob`` required two extra arguments when
it was called? By setting up our sampler with the ``args`` argument,
we're saying that the probability function should be called as:

.. code:: python
log_prob(p0[0], means, icov)
log_prob(p0[0], means, cov)
Expand All @@ -145,12 +138,6 @@ initial guess ``p0``:
pos, prob, state = sampler.run_mcmc(p0, 100)
sampler.reset()
.. parsed-literal::
100%|██████████| 100/100 [00:00<00:00, 305.62it/s]
You'll notice that I saved the final position of the walkers (after the
100 steps) to a variable called ``pos``. You can check out what will be
contained in the other output variables by looking at the documentation
Expand All @@ -166,56 +153,27 @@ Now, we can do our production run of 10000 steps:
sampler.run_mcmc(pos, 10000);
.. parsed-literal::
100%|██████████| 10000/10000 [00:24<00:00, 407.28it/s]
The sampler now has a property :attr:`EnsembleSampler.chain` that is a
numpy array with the shape ``(1000, 250, 50)``. Take note of that shape
and make sure that you know where each of those numbers come from.
Another useful object is the :attr:`EnsembleSampler.flatchain` which
has the shape ``(250000, 50)`` and contains all the samples reshaped
into a flat list. You can see now that we now have 250 000 unbiased
samples of the density :math:`p(\vec{x})`. You can make histograms of
these samples to get an estimate of the density that you were sampling:

:ref:`autocorr`

.. code:: python
sampler.get_autocorr_time()
.. parsed-literal::
array([ 54.61620237, 53.72668829, 54.67029465, 54.80001017, 53.99357549])
The samples can be accessed using the
:func:`EnsembleSampler.get_chain` method. This will return an array
with the shape ``(10000, 32, 5)`` giving the parameter values for each
walker at each step in the chain. Take note of that shape and make sure
that you know where each of those numbers come from. You can make
histograms of these samples to get an estimate of the density that you
were sampling:

.. code:: python
import matplotlib.pyplot as plt
for i in range(3):
plt.figure()
plt.hist(sampler.flatchain[:,i], 100, color="k", histtype="step")
plt.title("Dimension {0:d}".format(i))
.. image:: quickstart_files/quickstart_23_0.png


samples = sampler.get_chain(flat=True)
plt.hist(samples[:, 0], 100, color="k", histtype="step")
plt.xlabel(r"$\theta_1$")
plt.ylabel(r"$p(\theta_1)$")
plt.gca().set_yticks([]);
.. image:: quickstart_files/quickstart_23_1.png

.. image:: quickstart_files/quickstart_23_2.png
.. image:: quickstart_files/quickstart_21_0.png


Another good test of whether or not the sampling went well is to check
Expand All @@ -230,27 +188,20 @@ the mean acceptance fraction of the ensemble using the
.. parsed-literal::
Mean acceptance fraction: 0.551
Mean acceptance fraction: 0.553
This number should be between approximately 0.25 and 0.5 if everything
went as planned.
and the integrated autocorrelation time (see the :ref:`autocorr`
tutorial for more details)

.. code:: python
plt.plot(sampler.chain[:, :, 0]);
.. image:: quickstart_files/quickstart_27_0.png


.. code:: python
plt.plot(sampler.chain[:, :, -1]);
print("Mean autocorrelation time: {0:.3f} steps"
.format(np.mean(sampler.get_autocorr_time())))
.. parsed-literal::
.. image:: quickstart_files/quickstart_28_0.png
Mean autocorrelation time: 62.493 steps
49 changes: 3 additions & 46 deletions docs/user/faq.rst
Original file line number Diff line number Diff line change
Expand Up @@ -5,8 +5,6 @@ FAQ

**The not-so-frequently asked questions that still have useful answers**

.. _walkers:

What are "walkers"?
-------------------

Expand All @@ -24,55 +22,14 @@ The best technique seems to be to start in a small ball around the a priori
preferred position. Don't worry, the walkers quickly branch out and explore
the rest of the space.

Wrapping C++ code
-----------------

There are numerous ways to do it, see
`the python wiki
<https://wiki.python.org/moin/IntegratingPythonWithOtherLanguages#C.2FC.2B-.2B->`_.

Extra care has to be taken if mpi support is needed as the mpi4py module used by
emcee depends on the pickle module to send a function call to different
processors/cores.

A minimal extension of the mpi.py example in which the target density is coded
in C++ and wrapped with the `swig library <http://swig.org/>`_ is shown in this
`gist <https://gist.github.com/fredRos/7122649>`_. It also demonstrates the hacks
needed to get the pickling to work.


Parameter limits
----------------

In order to confine the walkers to a finite volume of the parameter space, have
your function return negative infinity outside of the volume corresponding to
the logarithm of 0 prior probability using::

return -numpy.inf

Note: if your function is written in C++, use::

return -std::numeric_limits<double>::infinity();

and avoid::

return -std::numeric_limits<double>::max();

as it does not have the desired effect.

Troubleshooting
---------------

**I'm getting weird spikes in my data/I have low acceptance fractions/both...
what should I do?**

Double the number of walkers. If that doesn't work, double it again. And
again. Until you run out of RAM. At that point, I don't know!

the logarithm of 0 prior probability using

**The walkers are getting stuck in "islands" of low likelihood. How can I
fix that?**
.. code-block:: python
Try increasing the number of walkers. If that doesn't work, you can try
pruning using a clustering algorithm like the one found in
`arxiv:1104.2612 <http://arxiv.org/abs/1104.2612>`_.
return -numpy.inf
4 changes: 0 additions & 4 deletions tests/unit/test_sampler.py
Original file line number Diff line number Diff line change
Expand Up @@ -26,10 +26,6 @@ def normal_log_prob(params):
])
)
def test_shapes(backend, moves, nwalkers=32, ndim=3, nsteps=100, seed=1234):
import emcee
print(emcee.__file__)
assert 0

# Set up the random number generator.
np.random.seed(seed)

Expand Down

0 comments on commit ff4e5ec

Please sign in to comment.