diff --git a/docs/conf.py b/docs/conf.py
index 169414e6c..4083653fa 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -1,3 +1,4 @@
+
 #!/usr/bin/env python3
 # -*- coding: utf-8 -*-
 #
@@ -28,7 +29,6 @@
 import projectq.setups.default
 import projectq.setups.grid
 import projectq.setups.ibm
-import projectq.setups.ibm16
 import projectq.setups.linear
 import projectq.setups.restrictedgateset
 import projectq.setups.decompositions
diff --git a/examples/bellpair_circuit.py b/examples/bellpair_circuit.py
index 108e224fa..96c001ffe 100755
--- a/examples/bellpair_circuit.py
+++ b/examples/bellpair_circuit.py
@@ -1,13 +1,20 @@
+import matplotlib.pyplot as plt
+
 from projectq import MainEngine
 from projectq.backends import CircuitDrawer
+from projectq.setups.default import get_engine_list
+from projectq.libs.hist import histogram
 
 from teleport import create_bell_pair
 
 # create a main compiler engine
 drawing_engine = CircuitDrawer()
-eng = MainEngine(drawing_engine)
+eng = MainEngine(engine_list = get_engine_list() + [drawing_engine])
 
-create_bell_pair(eng)
+qb0, qb1 = create_bell_pair(eng)
 
 eng.flush()
 print(drawing_engine.get_latex())
+
+histogram(eng.backend, [qb0, qb1])
+plt.show()
diff --git a/examples/ibm.py b/examples/ibm.py
index 11a81a832..eafd50b80 100755
--- a/examples/ibm.py
+++ b/examples/ibm.py
@@ -1,8 +1,11 @@
-import projectq.setups.ibm
+import matplotlib.pyplot as plt
+import getpass
+
+from projectq import MainEngine
 from projectq.backends import IBMBackend
+from projectq.libs.hist import histogram
 from projectq.ops import Measure, Entangle, All
-from projectq import MainEngine
-import getpass
+import projectq.setups.ibm
 
 
 def run_entangle(eng, num_qubits=3):
@@ -29,9 +32,12 @@ def run_entangle(eng, num_qubits=3):
     eng.flush()
 
     # access the probabilities via the back-end:
-    results = eng.backend.get_probabilities(qureg)
-    for state in results:
-        print("Measured {} with p = {}.".format(state, results[state]))
+    # results = eng.backend.get_probabilities(qureg)
+    # for state in results:
+    #     print("Measured {} with p = {}.".format(state, results[state]))
+    # or plot them directly:
+    histogram(eng.backend, qureg)
+    plt.show()
 
     # return one (random) measurement outcome.
     return [int(q) for q in qureg]
@@ -39,11 +45,11 @@ def run_entangle(eng, num_qubits=3):
 
 if __name__ == "__main__":
     #devices commonly available :
-    #ibmq_16_melbourne (15 qubit)
-    #ibmq_essex (5 qubit)
-    #ibmq_qasm_simulator (32 qubits)
-    device = None #replace by the IBM device name you want to use
-    token = None  #replace by the token given by IBMQ
+    # ibmq_16_melbourne (15 qubit)
+    # ibmq_essex (5 qubit)
+    # ibmq_qasm_simulator (32 qubits)
+    device = None # replace by the IBM device name you want to use
+    token = None  # replace by the token given by IBMQ
     if token is None:
         token = getpass.getpass(prompt='IBM Q token > ')
     if device is None:
diff --git a/examples/ibmq_tutorial.ipynb b/examples/ibmq_tutorial.ipynb
index 464612ca5..358bdcf9e 100644
--- a/examples/ibmq_tutorial.ipynb
+++ b/examples/ibmq_tutorial.ipynb
@@ -2,14 +2,18 @@
  "cells": [
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "# Running ProjectQ code on IBM Q devices"
    ]
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "In this tutorial, we will see how to run code on IBM Q devices directly from within ProjectQ. All that is needed is an IBM Q Experience user account. To sign up, visit https://quantumexperience.ng.bluemix.net/.\n",
     "\n",
@@ -19,10 +23,24 @@
     "First, we import all necessary operations (`Entangle`, measurement), the back-end (`IBMBackend`), and the main compiler engine (`MainEngine`). The Entangle operation is defined as a Hadamard gate on the first qubit (creates an equal superposition of |0> and |1>), followed by controlled NOT gates acting on all other qubits controlled on the first."
    ]
   },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {
+    "collapsed": false
+   },
+   "outputs": [],
+   "source": [
+    "%matplotlib inline\n",
+    "import matplotlib.pyplot as plt"
+   ]
+  },
   {
    "cell_type": "code",
    "execution_count": 1,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "import projectq.setups.ibm\n",
@@ -33,7 +51,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "Next, we instantiate a main compiler engine using the IBM Q back-end and the predefined compiler engines which take care of the qubit placement, translation of operations, etc.:"
    ]
@@ -41,7 +61,9 @@
   {
    "cell_type": "code",
    "execution_count": 2,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "eng = MainEngine(IBMBackend(use_hardware=True, num_runs=1024,\n",
@@ -51,7 +73,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "If `use_hardware` is set to `False`, it will use the IBM Q simulator instead. `num_runs` specifies the number of samples to collect for statistics, `verbose=True` would output additional information which may be helpful for debugging, and the device parameter lets users choose between the two devices (\"ibmqx4\" and \"ibmqx5\").\n",
     "\n",
@@ -61,7 +85,9 @@
   {
    "cell_type": "code",
    "execution_count": 3,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -105,9 +131,12 @@
     "    eng.flush()\n",
     "\n",
     "    # access the probabilities via the back-end:\n",
-    "    results = eng.backend.get_probabilities(qureg)\n",
-    "    for state in results:\n",
-    "        print(\"Measured {} with p = {}.\".format(state, results[state]))\n",
+    "    # results = eng.backend.get_probabilities(qureg)\n",
+    "    # for state in results:\n",
+    "    #     print(\"Measured {} with p = {}.\".format(state, results[state]))\n",
+    "    # or plot them directly:\n",
+    "    histogram(eng.backend, qureg)\n",
+    "    plt.show()\n",
     "\n",
     "    # return one (random) measurement outcome.\n",
     "    return [int(q) for q in qureg]\n",
@@ -117,7 +146,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "## Retrieving a timed-out execution\n",
     "Sometimes, the queue is very long and the waiting times may exceed the limit of 5 minutes. In this case, ProjectQ will raise an exception which contains the job ID, as could be seen above, where the job ID was `5b557df2306393003b746da2`."
@@ -125,7 +156,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "In order to still retrieve all results at a later point in time, one can simply re-run the entire program using a slightly modified back-end:"
    ]
@@ -133,7 +166,9 @@
   {
    "cell_type": "code",
    "execution_count": 4,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -193,7 +228,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "## Entangling more qubits: Using ibmqx5\n",
     "\n",
@@ -209,7 +246,9 @@
   {
    "cell_type": "code",
    "execution_count": 5,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [],
    "source": [
     "import projectq.setups.ibm16  # import setup which contains the grid mapper\n",
@@ -220,7 +259,9 @@
   },
   {
    "cell_type": "markdown",
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "source": [
     "and then re-run the example from before via `run_entangle(eng, num_qubits)`. If an execution times out, it can also be retrieved at a later point by providing the additional `retrieve_execution=\"execution_id\"` parameter to the IBMBackend (but this time with `device='ibmqx5'`)."
    ]
@@ -228,7 +269,9 @@
   {
    "cell_type": "code",
    "execution_count": 6,
-   "metadata": {},
+   "metadata": {
+    "collapsed": false
+   },
    "outputs": [
     {
      "name": "stdout",
@@ -488,8 +531,18 @@
  ],
  "metadata": {
   "kernelspec": {
+   "argv": [
+    "python",
+    "-m",
+    "ipykernel_launcher",
+    "-f",
+    "{connection_file}"
+   ],
    "display_name": "Python 3",
+   "env": null,
+   "interrupt_mode": "signal",
    "language": "python",
+   "metadata": null,
    "name": "python3"
   },
   "language_info": {
@@ -503,7 +556,8 @@
    "nbconvert_exporter": "python",
    "pygments_lexer": "ipython3",
    "version": "3.6.5"
-  }
+  },
+  "name": "ibmq_tutorial.ipynb"
  },
  "nbformat": 4,
  "nbformat_minor": 2
diff --git a/projectq/libs/hist/__init__.py b/projectq/libs/hist/__init__.py
new file mode 100644
index 000000000..088766263
--- /dev/null
+++ b/projectq/libs/hist/__init__.py
@@ -0,0 +1,20 @@
+#   Copyright 2020 ProjectQ-Framework (www.projectq.ch)
+#
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+"""
+contains a function to plot measurement outcome probabilities
+as a histogram for the simulator
+"""
+
+from ._histogram import histogram
diff --git a/projectq/libs/hist/_histogram.py b/projectq/libs/hist/_histogram.py
new file mode 100644
index 000000000..77b0d2c20
--- /dev/null
+++ b/projectq/libs/hist/_histogram.py
@@ -0,0 +1,78 @@
+#   Copyright 2020 ProjectQ-Framework (www.projectq.ch)
+#
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+from __future__ import print_function
+import matplotlib.pyplot as plt
+
+from projectq.backends import Simulator
+
+
+def histogram(backend, qureg):
+    """
+    Make a measurement outcome probability histogram for the given qubits.
+
+    Args:
+        backend (BasicEngine): A ProjectQ backend
+        qureg (list of qubits and/or quregs): The qubits,
+            for which to make the histogram
+
+    Returns:
+        A tuple (fig, axes, probabilities), where:
+        fig: The histogram as figure
+        axes: The axes of the histogram
+        probabilities (dict): A dictionary mapping outcomes as string
+            to their probabilities
+
+    Note:
+        Don't forget to call eng.flush() before using this function.
+    """
+    qubit_list = []
+    for q in qureg:
+        if isinstance(q, list):
+            qubit_list.extend(q)
+        else:
+            qubit_list.append(q)
+
+    if len(qubit_list) > 5:
+        print('Warning: For {0} qubits there are 2^{0} different outcomes'.
+              format(len(qubit_list)))
+        print("The resulting histogram may look bad and/or take too long.")
+        print("Consider calling histogram() with a sublist of the qubits.")
+
+    if hasattr(backend, 'get_probabilities'):
+        probabilities = backend.get_probabilities(qureg)
+    elif isinstance(backend, Simulator):
+        outcome = [0] * len(qubit_list)
+        n_outcomes = (1 << len(qubit_list))
+        probabilities = {}
+        for i in range(n_outcomes):
+            for pos in range(len(qubit_list)):
+                if (1 << pos) & i:
+                    outcome[pos] = 1
+                else:
+                    outcome[pos] = 0
+            probabilities[''.join([str(bit) for bit in outcome
+                                   ])] = backend.get_probability(
+                                       outcome, qubit_list)
+    else:
+        raise RuntimeError('Unable to retrieve probabilities from backend')
+
+    # Empirical figure size for up to 5 qubits
+    fig, axes = plt.subplots(figsize=(min(21.2, 2
+                                          + 0.6 * (1 << len(qubit_list))), 7))
+    names = list(probabilities.keys())
+    values = list(probabilities.values())
+    axes.bar(names, values)
+    fig.suptitle('Measurement Probabilities')
+    return (fig, axes, probabilities)
diff --git a/projectq/libs/hist/_histogram_test.py b/projectq/libs/hist/_histogram_test.py
new file mode 100644
index 000000000..1f38573b0
--- /dev/null
+++ b/projectq/libs/hist/_histogram_test.py
@@ -0,0 +1,135 @@
+#   Copyright 2020 ProjectQ-Framework (www.projectq.ch)
+#
+#   Licensed under the Apache License, Version 2.0 (the "License");
+#   you may not use this file except in compliance with the License.
+#   You may obtain a copy of the License at
+#
+#       http://www.apache.org/licenses/LICENSE-2.0
+#
+#   Unless required by applicable law or agreed to in writing, software
+#   distributed under the License is distributed on an "AS IS" BASIS,
+#   WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
+#   See the License for the specific language governing permissions and
+#   limitations under the License.
+
+import pytest
+import matplotlib
+import matplotlib.pyplot as plt
+
+from projectq import MainEngine
+from projectq.ops import H, C, X, Measure, All, AllocateQubitGate, FlushGate
+from projectq.cengines import DummyEngine, BasicEngine
+from projectq.backends import Simulator
+from projectq.libs.hist import histogram
+
+
+@pytest.fixture(scope="module")
+def matplotlib_setup():
+    old_backend = matplotlib.get_backend()
+    matplotlib.use('agg')  # avoid showing the histogram plots
+    yield
+    matplotlib.use(old_backend)
+
+
+def test_invalid_backend(matplotlib_setup):
+    eng = MainEngine(backend=DummyEngine())
+    qubit = eng.allocate_qubit()
+    eng.flush()
+
+    with pytest.raises(RuntimeError):
+        histogram(eng.backend, qubit)
+
+
+def test_backend_get_probabilities_method(matplotlib_setup):
+    class MyBackend(BasicEngine):
+        def get_probabilities(self, qureg):
+            return {'000': 0.5, '111': 0.5}
+
+        def is_available(self, cmd):
+            return True
+
+        def receive(self, command_list):
+            for cmd in command_list:
+                if not isinstance(cmd.gate, FlushGate):
+                    assert isinstance(cmd.gate, AllocateQubitGate)
+
+    eng = MainEngine(backend=MyBackend(), verbose=True)
+    qureg = eng.allocate_qureg(3)
+    eng.flush()
+    _, _, prob = histogram(eng.backend, qureg)
+    assert prob['000'] == 0.5
+    assert prob['111'] == 0.5
+
+
+def test_qubit(matplotlib_setup):
+    sim = Simulator()
+    eng = MainEngine(sim)
+    qubit = eng.allocate_qubit()
+    eng.flush()
+    _, _, prob = histogram(sim, qubit)
+    assert prob["0"] == pytest.approx(1)
+    assert prob["1"] == pytest.approx(0)
+    H | qubit
+    eng.flush()
+    _, _, prob = histogram(sim, qubit)
+    assert prob["0"] == pytest.approx(0.5)
+    Measure | qubit
+    eng.flush()
+    _, _, prob = histogram(sim, qubit)
+    assert prob["0"] == pytest.approx(1) or prob["1"] == pytest.approx(1)
+
+
+def test_qureg(matplotlib_setup):
+    sim = Simulator()
+    eng = MainEngine(sim)
+    qureg = eng.allocate_qureg(3)
+    eng.flush()
+    _, _, prob = histogram(sim, qureg)
+    assert prob["000"] == pytest.approx(1)
+    assert prob["110"] == pytest.approx(0)
+    H | qureg[0]
+    C(X, 1) | (qureg[0], qureg[1])
+    H | qureg[2]
+    eng.flush()
+    _, _, prob = histogram(sim, qureg)
+    assert prob["110"] == pytest.approx(0.25)
+    assert prob["100"] == pytest.approx(0)
+    All(Measure) | qureg
+    eng.flush()
+    _, _, prob = histogram(sim, qureg)
+    assert prob["000"] == pytest.approx(1) or prob["001"] == pytest.approx(1) \
+        or prob["110"] == pytest.approx(1) or prob["111"] == pytest.approx(1)
+    assert prob["000"] + prob["001"] + prob["110"] + prob[
+        "111"] == pytest.approx(1)
+
+
+def test_combination(matplotlib_setup):
+    sim = Simulator()
+    eng = MainEngine(sim)
+    qureg = eng.allocate_qureg(2)
+    qubit = eng.allocate_qubit()
+    eng.flush()
+    _, _, prob = histogram(sim, [qureg, qubit])
+    assert prob["000"] == pytest.approx(1)
+    H | qureg[0]
+    C(X, 1) | (qureg[0], qureg[1])
+    H | qubit
+    Measure | qureg[0]
+    eng.flush()
+    _, _, prob = histogram(sim, [qureg, qubit])
+    assert (prob["000"] == pytest.approx(0.5) and prob["001"] == pytest.approx(0.5)) \
+        or (prob["110"] == pytest.approx(0.5) and prob["111"] == pytest.approx(0.5))
+    assert prob["100"] == pytest.approx(0)
+    Measure | qubit
+
+
+def test_too_many_qubits(matplotlib_setup, capsys):
+    sim = Simulator()
+    eng = MainEngine(sim)
+    qureg = eng.allocate_qureg(6)
+    eng.flush()
+    l_ref = len(capsys.readouterr().out)
+    _, _, prob = histogram(sim, qureg)
+    assert len(capsys.readouterr().out) > l_ref
+    assert prob["000000"] == pytest.approx(1)
+    All(Measure)
diff --git a/pytest.ini b/pytest.ini
index fab634b12..d17d4ce2e 100755
--- a/pytest.ini
+++ b/pytest.ini
@@ -4,3 +4,4 @@ testpaths = projectq
 filterwarnings =
     error
     ignore:the matrix subclass is not the recommended way:PendingDeprecationWarning
+    ignore:Matplotlib is currently using agg, which is a non-GUI backend, so cannot show the figure.