diff --git a/docs/conf.py b/docs/conf.py
index 4eee9a879..5516fbf3e 100755
--- a/docs/conf.py
+++ b/docs/conf.py
@@ -119,8 +119,8 @@
     "examples/compare_tempo2_B1855.py",
     "examples/example_dmx_ranges.py",
     "examples/example_pulse_numbers.py",
-    "examples/bayesian-example-NGC6440E.py",
-    "examples/bayesian-wideband-example.py",
+    # "examples/bayesian-example-NGC6440E.py",
+    # "examples/bayesian-wideband-example.py",
     "conf.py",
     "_ext",
 ]
diff --git a/docs/examples/bayesian-example-NGC6440E.py b/docs/examples/bayesian-example-NGC6440E.py
index f19afd219..d371c7316 100644
--- a/docs/examples/bayesian-example-NGC6440E.py
+++ b/docs/examples/bayesian-example-NGC6440E.py
@@ -75,31 +75,41 @@
     + np.random.randn(nwalkers, bt.nparams) * maxlike_errors
 )
 
+# %%
+# ** IMPORTANT!!! **
+# This is used to exclude some of the following time-consuming steps from the readthedocs build.
+# Set this to False while actually using this example.
+rtd = True
+
 # %%
 # Use longer chain_length for real runs. It is kept small here so that
 # the sampling finishes quickly (and because I know the burn in is short
 # because of the cheating priors above).
-print("Running emcee...")
-chain_length = 1000
-sampler.run_mcmc(
-    start_points,
-    chain_length,
-    progress=True,
-)
+if not rtd:
+    print("Running emcee...")
 
-# %%
-# Merge all the chains together after discarding the first 100 samples as 'burn-in'.
-# The burn-in should be decided after looking at the chains in the real world.
-samples_emcee = sampler.get_chain(flat=True, discard=100)
+    chain_length = 1000
+
+    sampler.run_mcmc(
+        start_points,
+        chain_length,
+        progress=True,
+    )
 
 # %%
-# Plot the MCMC chains to make sure that the burn-in has been removed properly.
-# Otherwise, go back and discard more points.
-for idx, param_chain in enumerate(samples_emcee.T):
-    plt.subplot(bt.nparams, 1, idx + 1)
-    plt.plot(param_chain, label=bt.param_labels[idx])
-    plt.legend()
-plt.show()
+if not rtd:
+    # Merge all the chains together after discarding the first 100 samples as 'burn-in'.
+    # The burn-in should be decided after looking at the chains in the real world.
+    samples_emcee = sampler.get_chain(flat=True, discard=100)
+
+    # %%
+    # Plot the MCMC chains to make sure that the burn-in has been removed properly.
+    # Otherwise, go back and discard more points.
+    for idx, param_chain in enumerate(samples_emcee.T):
+        plt.subplot(bt.nparams, 1, idx + 1)
+        plt.plot(param_chain, label=bt.param_labels[idx])
+        plt.legend()
+    plt.show()
 
 # %%
 # Plot the posterior distribution.
@@ -123,28 +133,30 @@
 # `dlogz` is the target accuracy in the computed Bayesian evidence.
 # Increasing `npoints` or decreasing `dlogz` gives more accurate results,
 # but at the cost of time.
-print("Running nestle...")
-result_nestle_1 = nestle.sample(
-    bt.lnlikelihood,
-    bt.prior_transform,
-    bt.nparams,
-    method="multi",
-    npoints=150,
-    dlogz=0.5,
-    callback=nestle.print_progress,
-)
+if not rtd:
+    print("Running nestle...")
+    result_nestle_1 = nestle.sample(
+        bt.lnlikelihood,
+        bt.prior_transform,
+        bt.nparams,
+        method="multi",
+        npoints=150,
+        dlogz=0.5,
+        callback=nestle.print_progress,
+    )
 
 # %%
 # Plot the posterior
 # The nested samples come with weights, which must be taken into account
 # while plotting.
-fig = corner.corner(
-    result_nestle_1.samples,
-    weights=result_nestle_1.weights,
-    labels=bt.param_labels,
-    range=[0.999] * bt.nparams,
-)
-plt.show()
+if not rtd:
+    fig = corner.corner(
+        result_nestle_1.samples,
+        weights=result_nestle_1.weights,
+        labels=bt.param_labels,
+        range=[0.999] * bt.nparams,
+    )
+    plt.show()
 
 # %% [markdown]
 # Let us create a new model with an EFAC applied to all toas (all
@@ -175,36 +187,41 @@
 print(bt2.likelihood_method)
 
 # %%
-result_nestle_2 = nestle.sample(
-    bt2.lnlikelihood,
-    bt2.prior_transform,
-    bt2.nparams,
-    method="multi",
-    npoints=150,
-    dlogz=0.5,
-    callback=nestle.print_progress,
-)
+if not rtd:
+    result_nestle_2 = nestle.sample(
+        bt2.lnlikelihood,
+        bt2.prior_transform,
+        bt2.nparams,
+        method="multi",
+        npoints=150,
+        dlogz=0.5,
+        callback=nestle.print_progress,
+    )
 
 # %%
 # Plot the posterior.
 # The EFAC looks consistent with 1.
-fig2 = corner.corner(
-    result_nestle_2.samples,
-    weights=result_nestle_2.weights,
-    labels=bt2.param_labels,
-    range=[0.999] * bt2.nparams,
-)
-plt.show()
+if not rtd:
+    fig2 = corner.corner(
+        result_nestle_2.samples,
+        weights=result_nestle_2.weights,
+        labels=bt2.param_labels,
+        range=[0.999] * bt2.nparams,
+    )
+    plt.show()
 
 # %% [markdown]
 # Now let us look at the evidences and compute the Bayes factor.
 
 # %%
-print(f"Evidence without EFAC : {result_nestle_1.logz} +/- {result_nestle_1.logzerr}")
-print(f"Evidence with EFAC : {result_nestle_2.logz} +/- {result_nestle_2.logzerr}")
+if not rtd:
+    print(
+        f"Evidence without EFAC : {result_nestle_1.logz} +/- {result_nestle_1.logzerr}"
+    )
+    print(f"Evidence with EFAC : {result_nestle_2.logz} +/- {result_nestle_2.logzerr}")
 
-bf = np.exp(result_nestle_1.logz - result_nestle_2.logz)
-print(f"Bayes factor : {bf} (in favor of no EFAC)")
+    bf = np.exp(result_nestle_1.logz - result_nestle_2.logz)
+    print(f"Bayes factor : {bf} (in favor of no EFAC)")
 
 # %% [markdown]
 # The Bayes factor tells us that the EFAC is unnecessary for this dataset.
diff --git a/docs/examples/bayesian-wideband-example.py b/docs/examples/bayesian-wideband-example.py
index 85f949eea..f576b9546 100644
--- a/docs/examples/bayesian-wideband-example.py
+++ b/docs/examples/bayesian-wideband-example.py
@@ -72,28 +72,36 @@
 )
 
 # %%
-print("Running emcee...")
-chain_length = 1000
-sampler.run_mcmc(
-    start_points,
-    chain_length,
-    progress=True,
-)
+# ** IMPORTANT!!! **
+# This is used to exclude the following time-consuming steps from the readthedocs build.
+# Set this to False while actually using this example.
+rtd = True
 
 # %%
-samples_emcee = sampler.get_chain(flat=True, discard=100)
+if not rtd:
+    print("Running emcee...")
+    chain_length = 1000
+    sampler.run_mcmc(
+        start_points,
+        chain_length,
+        progress=True,
+    )
+
+    samples_emcee = sampler.get_chain(flat=True, discard=100)
 
 # %%
 # Plot the chains to make sure they have converged and the burn-in has been removed properly.
-for idx, param_chain in enumerate(samples_emcee.T):
-    plt.subplot(bt.nparams, 1, idx + 1)
-    plt.plot(param_chain)
-    plt.ylabel(bt.param_labels[idx])
-    plt.autoscale()
-plt.show()
+if not rtd:
+    for idx, param_chain in enumerate(samples_emcee.T):
+        plt.subplot(bt.nparams, 1, idx + 1)
+        plt.plot(param_chain)
+        plt.ylabel(bt.param_labels[idx])
+        plt.autoscale()
+    plt.show()
 
 # %%
-fig = corner.corner(
-    samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params
-)
-plt.show()
+if not rtd:
+    fig = corner.corner(
+        samples_emcee, labels=bt.param_labels, quantiles=[0.5], truths=maxlike_params
+    )
+    plt.show()