diff --git a/bilby/core/result.py b/bilby/core/result.py
index 99efa4912244595ffb2db0ab5e125234bea3eb59..48797a4f100e6583a880b59822df61997c6aeaa8 100644
--- a/bilby/core/result.py
+++ b/bilby/core/result.py
@@ -1994,9 +1994,9 @@ def plot_multiple(results, filename=None, labels=None, colours=None,
 
     if evidences:
         if np.isnan(results[0].log_bayes_factor):
-            template = ' $\mathrm{{ln}}(Z)={lnz:1.3g}$'
+            template = r' $\mathrm{{ln}}(Z)={lnz:1.3g}$'
         else:
-            template = ' $\mathrm{{ln}}(B)={lnbf:1.3g}$'
+            template = r' $\mathrm{{ln}}(B)={lnbf:1.3g}$'
         labels = [template.format(lnz=result.log_evidence,
                                   lnbf=result.log_bayes_factor)
                   for ii, result in enumerate(results)]
diff --git a/bilby/gw/likelihood/base.py b/bilby/gw/likelihood/base.py
index d6c9f6f14e0e74de4cdd20aec6e0e22175a3aba4..b468ff67ed793a9644288b08d5b0b0758d37ca49 100644
--- a/bilby/gw/likelihood/base.py
+++ b/bilby/gw/likelihood/base.py
@@ -651,7 +651,7 @@ class GravitationalWaveTransient(Likelihood):
 
     def generate_phase_sample_from_marginalized_likelihood(
             self, signal_polarizations=None):
-        """
+        r"""
         Generate a single sample from the posterior distribution for phase when
         using a likelihood which explicitly marginalises over phase.
 
diff --git a/bilby/gw/likelihood/multiband.py b/bilby/gw/likelihood/multiband.py
index 763912e84da1624206a233a8deb53357e59c8190..f361742c123af9cd5f27d236b76d1f3e1a4efc7b 100644
--- a/bilby/gw/likelihood/multiband.py
+++ b/bilby/gw/likelihood/multiband.py
@@ -216,7 +216,7 @@ class MBGravitationalWaveTransient(GravitationalWaveTransient):
 
     @maximum_banding_frequency.setter
     def maximum_banding_frequency(self, maximum_banding_frequency):
-        """
+        r"""
         This sets the upper limit on a starting frequency of a band. The default value is the frequency at
         which f - 1 / \sqrt(- d\tau / df) starts to decrease, because the bisection search of the starting
         frequency does not work from that frequency. The stationary phase approximation is not valid at such
@@ -341,7 +341,7 @@ class MBGravitationalWaveTransient(GravitationalWaveTransient):
         return f, 1. / np.sqrt(-self._dtaudf(f))
 
     def _setup_frequency_bands(self):
-        """Set up frequency bands. The durations of bands geometrically decrease T, T/2. T/4, ..., where T is the
+        r"""Set up frequency bands. The durations of bands geometrically decrease T, T/2. T/4, ..., where T is the
         original duration. This sets the following instance variables.
 
         durations: durations of bands (T^(b) in the paper)
diff --git a/bilby/gw/prior.py b/bilby/gw/prior.py
index 547759740018bd3810e5d6baa50fba92c7045bf2..1bf3c69981d91b2edb10fa808b0ce5ffbeecae06 100644
--- a/bilby/gw/prior.py
+++ b/bilby/gw/prior.py
@@ -68,7 +68,7 @@ def convert_to_flat_in_component_mass_prior(result, fraction=0.25):
         priors[key] = Constraint(priors[key].minimum, priors[key].maximum, key, latex_label=priors[key].latex_label)
     for key in ['mass_1', 'mass_2']:
         priors[key] = Uniform(priors[key].minimum, priors[key].maximum, key, latex_label=priors[key].latex_label,
-                              unit="$M_{\odot}$")
+                              unit=r"$M_{\odot}$")
 
     weights = np.array(result.get_weights_by_new_prior(old_priors, priors,
                                                        prior_names=['chirp_mass', 'mass_ratio', 'mass_1', 'mass_2']))
@@ -276,7 +276,7 @@ class UniformComovingVolume(Cosmological):
 
 
 class UniformSourceFrame(Cosmological):
-    """
+    r"""
     Prior for redshift which is uniform in comoving volume and source frame
     time.
 
@@ -294,7 +294,7 @@ class UniformSourceFrame(Cosmological):
 class UniformInComponentsChirpMass(PowerLaw):
 
     def __init__(self, minimum, maximum, name='chirp_mass',
-                 latex_label='$\mathcal{M}$', unit=None, boundary=None):
+                 latex_label=r'$\mathcal{M}$', unit=None, boundary=None):
         """
         Prior distribution for chirp mass which is uniform in component masses.
 
@@ -924,36 +924,36 @@ Prior._default_latex_labels = {
     'mass_1': '$m_1$',
     'mass_2': '$m_2$',
     'total_mass': '$M$',
-    'chirp_mass': '$\mathcal{M}$',
+    'chirp_mass': r'$\mathcal{M}$',
     'mass_ratio': '$q$',
-    'symmetric_mass_ratio': '$\eta$',
+    'symmetric_mass_ratio': r'$\eta$',
     'a_1': '$a_1$',
     'a_2': '$a_2$',
-    'tilt_1': '$\\theta_1$',
-    'tilt_2': '$\\theta_2$',
-    'cos_tilt_1': '$\cos\\theta_1$',
-    'cos_tilt_2': '$\cos\\theta_2$',
-    'phi_12': '$\Delta\phi$',
-    'phi_jl': '$\phi_{JL}$',
+    'tilt_1': r'$\theta_1$',
+    'tilt_2': r'$\theta_2$',
+    'cos_tilt_1': r'$\cos\theta_1$',
+    'cos_tilt_2': r'$\cos\theta_2$',
+    'phi_12': r'$\Delta\phi$',
+    'phi_jl': r'$\phi_{JL}$',
     'luminosity_distance': '$d_L$',
-    'dec': '$\mathrm{DEC}$',
-    'ra': '$\mathrm{RA}$',
-    'iota': '$\iota$',
-    'cos_iota': '$\cos\iota$',
-    'theta_jn': '$\\theta_{JN}$',
-    'cos_theta_jn': '$\cos\\theta_{JN}$',
-    'psi': '$\psi$',
-    'phase': '$\phi$',
+    'dec': r'$\mathrm{DEC}$',
+    'ra': r'$\mathrm{RA}$',
+    'iota': r'$\iota$',
+    'cos_iota': r'$\cos\iota$',
+    'theta_jn': r'$\theta_{JN}$',
+    'cos_theta_jn': r'$\cos\theta_{JN}$',
+    'psi': r'$\psi$',
+    'phase': r'$\phi$',
     'geocent_time': '$t_c$',
     'time_jitter': '$t_j$',
-    'lambda_1': '$\\Lambda_1$',
-    'lambda_2': '$\\Lambda_2$',
-    'lambda_tilde': '$\\tilde{\\Lambda}$',
-    'delta_lambda_tilde': '$\\delta\\tilde{\\Lambda}$',
-    'chi_1': '$\\chi_1$',
-    'chi_2': '$\\chi_2$',
-    'chi_1_in_plane': '$\\chi_{1, \perp}$',
-    'chi_2_in_plane': '$\\chi_{2, \perp}$',
+    'lambda_1': r'$\Lambda_1$',
+    'lambda_2': r'$\Lambda_2$',
+    'lambda_tilde': r'$\tilde{\Lambda}$',
+    'delta_lambda_tilde': r'$\delta\tilde{\Lambda}$',
+    'chi_1': r'$\chi_1$',
+    'chi_2': r'$\chi_2$',
+    'chi_1_in_plane': r'$\chi_{1, \perp}$',
+    'chi_2_in_plane': r'$\chi_{2, \perp}$',
 }
 
 
@@ -1056,7 +1056,7 @@ class CalibrationPriorDict(PriorDict):
                                    boundary=boundary)
         for ii in range(n_nodes):
             name = "recalib_{}_phase_{}".format(label, ii)
-            latex_label = "$\\phi^{}_{}$".format(label, ii)
+            latex_label = r"$\phi^{}_{}$".format(label, ii)
             prior[name] = Gaussian(mu=phase_mean_nodes[ii],
                                    sigma=phase_sigma_nodes[ii],
                                    name=name, latex_label=latex_label,
@@ -1117,7 +1117,7 @@ class CalibrationPriorDict(PriorDict):
                                    boundary='reflective')
         for ii in range(n_nodes):
             name = "recalib_{}_phase_{}".format(label, ii)
-            latex_label = "$\\phi^{}_{}$".format(label, ii)
+            latex_label = r"$\phi^{}_{}$".format(label, ii)
             prior[name] = Gaussian(mu=phase_mean_nodes[ii],
                                    sigma=phase_sigma_nodes[ii],
                                    name=name, latex_label=latex_label,
diff --git a/bilby/gw/result.py b/bilby/gw/result.py
index fc599232b3a8cf71633825a4c941b61347b12ba3..5c95fe0f149625ca9360e6fbe5c8aa21d83c46cd 100644
--- a/bilby/gw/result.py
+++ b/bilby/gw/result.py
@@ -185,7 +185,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
             if len(amp_params) > 0:
                 amplitude = 100 * np.column_stack([posterior[param] for param in amp_params])
                 plot_spline_pos(logfreqs, amplitude, color=color, level=level,
-                                label="{0} (mean, {1}$\%$)".format(ifo.upper(), int(level * 100)))
+                                label=r"{0} (mean, {1}$\%$)".format(ifo.upper(), int(level * 100)))
 
             # Phase calibration model
             plt.sca(ax2)
@@ -194,7 +194,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
             if len(phase_params) > 0:
                 phase = np.column_stack([posterior[param] for param in phase_params])
                 plot_spline_pos(logfreqs, phase, color=color, level=level,
-                                label="{0} (mean, {1}$\%$)".format(ifo.upper(), int(level * 100)),
+                                label=r"{0} (mean, {1}$\%$)".format(ifo.upper(), int(level * 100)),
                                 xform=spline_angle_xform)
 
         ax1.tick_params(labelsize=.75 * font_size)
@@ -204,7 +204,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
         ax2.set_xscale('log')
 
         ax2.set_xlabel('Frequency [Hz]', fontsize=font_size)
-        ax1.set_ylabel('Amplitude [$\%$]', fontsize=font_size)
+        ax1.set_ylabel(r'Amplitude [$\%$]', fontsize=font_size)
         ax2.set_ylabel('Phase [deg]', fontsize=font_size)
 
         filename = os.path.join(outdir, self.label + '_calibration.' + format)
@@ -585,7 +585,7 @@ class CompactBinaryCoalescenceResult(CoreResult):
                 plot_frequencies,
                 np.percentile(fd_waveforms, lower_percentile, axis=0),
                 np.percentile(fd_waveforms, upper_percentile, axis=0),
-                color=WAVEFORM_COLOR, label='{}\% credible interval'.format(
+                color=WAVEFORM_COLOR, label=r'{}\% credible interval'.format(
                     int(upper_percentile - lower_percentile)),
                 alpha=0.3)
             axs[1].plot(
diff --git a/examples/core_examples/hyper_parameter_example.py b/examples/core_examples/hyper_parameter_example.py
index 0798a9faa07315d1f869d8bf8c6ff09fe0fc0e4a..d9ddba4b72b086f8c5d4b51f7967a1ccc78afe16 100644
--- a/examples/core_examples/hyper_parameter_example.py
+++ b/examples/core_examples/hyper_parameter_example.py
@@ -76,8 +76,8 @@ hp_likelihood = HyperparameterLikelihood(
     posteriors=samples, hyper_prior=hyper_prior,
     sampling_prior=run_prior, log_evidences=evidences, max_samples=500)
 
-hp_priors = dict(mu=Uniform(-10, 10, 'mu', '$\mu_{c0}$'),
-                 sigma=Uniform(0, 10, 'sigma', '$\sigma_{c0}$'))
+hp_priors = dict(mu=Uniform(-10, 10, 'mu', r'$\mu_{c0}$'),
+                 sigma=Uniform(0, 10, 'sigma', r'$\sigma_{c0}$'))
 
 # And run sampler
 result = run_sampler(
diff --git a/examples/core_examples/slabspike_example.py b/examples/core_examples/slabspike_example.py
index d798512353178bf9f77940c0017a2e1f21aa1b35..d7501580bf8dc880493d9f54908b17f187697180 100644
--- a/examples/core_examples/slabspike_example.py
+++ b/examples/core_examples/slabspike_example.py
@@ -64,12 +64,12 @@ priors['amplitude_2'] = bilby.core.prior.SlabSpikePrior(slab=amplitude_slab_2, s
 # Our problem has a degeneracy in the ordering. In general, this problem is somewhat difficult to resolve properly.
 # See e.g. https://github.com/GregoryAshton/kookaburra/blob/master/src/priors.py#L72 for an implementation.
 # We resolve this by not letting the priors overlap in this case.
-priors['mu_0'] = bilby.core.prior.Uniform(minimum=-5, maximum=-2, name='mu_0', latex_label='$\mu_0$')
-priors['mu_1'] = bilby.core.prior.Uniform(minimum=-2, maximum=2, name='mu_1', latex_label='$\mu_1$')
-priors['mu_2'] = bilby.core.prior.Uniform(minimum=2, maximum=5, name='mu_2', latex_label='$\mu_2$')
-priors['sigma_0'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_0', latex_label='$\sigma_0$')
-priors['sigma_1'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_1', latex_label='$\sigma_1$')
-priors['sigma_2'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_2', latex_label='$\sigma_2$')
+priors['mu_0'] = bilby.core.prior.Uniform(minimum=-5, maximum=-2, name='mu_0', latex_label=r'$\mu_0$')
+priors['mu_1'] = bilby.core.prior.Uniform(minimum=-2, maximum=2, name='mu_1', latex_label=r'$\mu_1$')
+priors['mu_2'] = bilby.core.prior.Uniform(minimum=2, maximum=5, name='mu_2', latex_label=r'$\mu_2$')
+priors['sigma_0'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_0', latex_label=r'$\sigma_0$')
+priors['sigma_1'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_1', latex_label=r'$\sigma_1$')
+priors['sigma_2'] = bilby.core.prior.LogUniform(minimum=0.01, maximum=10, name='sigma_2', latex_label=r'$\sigma_2$')
 
 # Setting up the likelihood and running the samplers works the same as elsewhere.
 likelihood = bilby.core.likelihood.GaussianLikelihood(x=xs, y=ys, func=triple_gaussian, sigma=sigma)
diff --git a/setup.cfg b/setup.cfg
index 4a86df8ebe1f77c62f4e6de88f45610732c0b3fe..7f999dcf7d5ee889b4a528f229d1d17801451a55 100644
--- a/setup.cfg
+++ b/setup.cfg
@@ -1,7 +1,7 @@
 [flake8]
 exclude = .git,docs,build,dist,test,*__init__.py
 max-line-length = 120
-ignore = E129 W503 W504 W605 E203 E402
+ignore = E129 W503 W504 E203 E402
 
 [tool:pytest]
 addopts =
diff --git a/test/core/prior/base_test.py b/test/core/prior/base_test.py
index bafa10fa90e7495ccf0999f900580aee65aa719a..49887de5032580096b1acdf321f5a83534ab5d3c 100644
--- a/test/core/prior/base_test.py
+++ b/test/core/prior/base_test.py
@@ -100,7 +100,7 @@ class TestPrior(unittest.TestCase):
     def test_default_label_assignment(self):
         self.prior.name = "chirp_mass"
         self.prior.latex_label = None
-        self.assertEqual(self.prior.latex_label, "$\mathcal{M}$")
+        self.assertEqual(self.prior.latex_label, r"$\mathcal{M}$")
 
     def test_default_label_assignment_default(self):
         self.assertTrue(self.prior.latex_label, self.prior.name)
diff --git a/test/core/prior/dict_test.py b/test/core/prior/dict_test.py
index 3b8487fa1f49ea9c25459c738cc16e58252a5f76..03881612a61404db58ebc43d2944fc833480643b 100644
--- a/test/core/prior/dict_test.py
+++ b/test/core/prior/dict_test.py
@@ -66,7 +66,7 @@ class TestPriorDict(unittest.TestCase):
                 name="chirp_mass",
                 minimum=25,
                 maximum=100,
-                latex_label="$\mathcal{M}$",
+                latex_label=r"$\mathcal{M}$",
             ),
             mass_ratio=bilby.core.prior.Uniform(
                 name="mass_ratio",
@@ -176,7 +176,7 @@ class TestPriorDict(unittest.TestCase):
                     name="chirp_mass",
                     minimum=25,
                     maximum=100,
-                    latex_label="$\mathcal{M}$",
+                    latex_label=r"$\mathcal{M}$",
                 ),
                 mass_ratio=bilby.core.prior.Uniform(
                     name="mass_ratio",