diff --git a/pyerrors/correlators.py b/pyerrors/correlators.py
index 10bc3b14..01a7d61e 100644
--- a/pyerrors/correlators.py
+++ b/pyerrors/correlators.py
@@ -1,3 +1,4 @@
+import warnings
 import numpy as np
 import autograd.numpy as anp
 import matplotlib.pyplot as plt
@@ -130,7 +131,7 @@ class Corr:
             raise Exception("Can not symmetrize odd T")
 
         if np.argmax(np.abs(self.content)) != 0:
-            print('Warning: correlator does not seem to be symmetric around x0=0.')
+            warnings.warn("Correlator does not seem to be symmetric around x0=0.", RuntimeWarning)
 
         newcontent = [self.content[0]]
         for t in range(1, self.T):
@@ -148,7 +149,7 @@ class Corr:
             raise Exception("Can not symmetrize odd T")
 
         if not all([o.zero_within_error() for o in self.content[0]]):
-            print('Warning: correlator does not seem to be anti-symmetric around x0=0.')
+            warnings.warn("Correlator does not seem to be anti-symmetric around x0=0.", RuntimeWarning)
 
         newcontent = [self.content[0]]
         for t in range(1, self.T):
diff --git a/pyerrors/fits.py b/pyerrors/fits.py
index 137ede2f..da8aed6b 100644
--- a/pyerrors/fits.py
+++ b/pyerrors/fits.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 # coding: utf-8
 
+import warnings
 import numpy as np
 import autograd.numpy as anp
 import scipy.optimize
@@ -300,7 +301,7 @@ def odr_fit(x, y, func, silent=False, **kwargs):
         P_phi = A @ np.linalg.inv(A.T @ A) @ A.T
         expected_chisquare = np.trace((np.identity(P_phi.shape[0]) - P_phi) @ W @ cov @ W)
         if expected_chisquare <= 0.0:
-            print('Warning, negative expected_chisquare.')
+            warnings.warn("Negative expected_chisquare.", RuntimeWarning)
             expected_chisquare = np.abs(expected_chisquare)
         result_dict['chisquare/expected_chisquare'] = odr_chisquare(np.concatenate((output.beta, output.xplus.ravel()))) / expected_chisquare
         if not silent:
@@ -611,9 +612,7 @@ def error_band(x, func, beta):
     """Returns the error band for an array of sample values x, for given fit function func with optimized parameters beta."""
     cov = covariance_matrix(beta)
     if np.any(np.abs(cov - cov.T) > 1000 * np.finfo(np.float64).eps):
-        print('Warning, Covariance matrix is not symmetric within floating point precision')
-        print('cov - cov.T:')
-        print(cov - cov.T)
+        warnings.warn("Covariance matrix is not symmetric within floating point precision", RuntimeWarning)
 
     deriv = []
     for i, item in enumerate(x):
diff --git a/pyerrors/pyerrors.py b/pyerrors/pyerrors.py
index 9a0d02c0..15d01f7f 100644
--- a/pyerrors/pyerrors.py
+++ b/pyerrors/pyerrors.py
@@ -1,6 +1,7 @@
 #!/usr/bin/env python
 # coding: utf-8
 
+import warnings
 import pickle
 import numpy as np
 import autograd.numpy as anp  # Thinly-wrapped numpy
@@ -784,9 +785,9 @@ def correlate(obs_a, obs_b):
             raise Exception('Shapes of ensemble', name, 'do not fit')
 
     if obs_a.reweighted == 1:
-        print('Warning: The first observable is already reweighted.')
+        warnings.warn("The first observable is already reweighted.", RuntimeWarning)
     if obs_b.reweighted == 1:
-        print('Warning: The second observable is already reweighted.')
+        warnings.warn("The second observable is already reweighted.", RuntimeWarning)
 
     new_samples = []
     for name in sorted(obs_a.names):