diff --git a/GPy/examples/coreg_example.py b/GPy/examples/coreg_example.py
index 4e9566dc..3afeb9fb 100644
--- a/GPy/examples/coreg_example.py
+++ b/GPy/examples/coreg_example.py
@@ -3,7 +3,7 @@
 
 import numpy as np
 try:
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 import GPy
diff --git a/GPy/examples/non_gaussian.py b/GPy/examples/non_gaussian.py
index 3652b4d3..5d5fc7f6 100644
--- a/GPy/examples/non_gaussian.py
+++ b/GPy/examples/non_gaussian.py
@@ -77,7 +77,7 @@ def student_t_approx(optimize=True, plot=True):
     debug=True
     if debug:
         m4.optimize(messages=1)
-        import pylab as pb
+        from matplotlib import pyplot as pb
         pb.plot(m4.X, m4.inference_method.f_hat)
         pb.plot(m4.X, m4.Y, 'rx')
         m4.plot()
diff --git a/GPy/examples/regression.py b/GPy/examples/regression.py
index 267c6d1e..1ce2f6f4 100644
--- a/GPy/examples/regression.py
+++ b/GPy/examples/regression.py
@@ -5,7 +5,7 @@
 Gaussian Processes regression examples
 """
 try:
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 import numpy as np
diff --git a/GPy/models/gradient_checker.py b/GPy/models/gradient_checker.py
index 728f3e23..42d41e60 100644
--- a/GPy/models/gradient_checker.py
+++ b/GPy/models/gradient_checker.py
@@ -251,7 +251,7 @@ class HessianChecker(GradientChecker):
             print(grad_string)
 
             if plot:
-                import pylab as pb
+                from matplotlib import pyplot as pb
                 fig, axes = pb.subplots(2, 2)
                 max_lim = numpy.max(numpy.vstack((analytic_hess, numeric_hess)))
                 min_lim = numpy.min(numpy.vstack((analytic_hess, numeric_hess)))
diff --git a/GPy/plotting/__init__.py b/GPy/plotting/__init__.py
index 9dd84441..dc15724e 100644
--- a/GPy/plotting/__init__.py
+++ b/GPy/plotting/__init__.py
@@ -4,4 +4,4 @@
 try:
     from . import matplot_dep
 except (ImportError, NameError):
-    print('Fail to load GPy.plotting.matplot_dep.')
\ No newline at end of file
+    # Matplotlib not available
\ No newline at end of file
diff --git a/GPy/plotting/matplot_dep/Tango.py b/GPy/plotting/matplot_dep/Tango.py
index eeb2e075..5c004519 100644
--- a/GPy/plotting/matplot_dep/Tango.py
+++ b/GPy/plotting/matplot_dep/Tango.py
@@ -3,7 +3,7 @@
 
 
 import matplotlib as mpl
-import pylab as pb
+from matplotlib import pyplot as pb
 import sys
 #sys.path.append('/home/james/mlprojects/sitran_cluster/')
 #from switch_pylab_backend import *
@@ -159,7 +159,7 @@ cdict_Alu = {'red' :((0./5,colorsRGB['Aluminium1'][0]/256.,colorsRGB['Aluminium1
 # cmap_BGR = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_BGR,256)
 # cmap_RB = mpl.colors.LinearSegmentedColormap('TangoRedBlue',cdict_RB,256)
 if __name__=='__main__':
-    import pylab as pb
+    from matplotlib import pyplot as pb
     pb.figure()
     pb.pcolor(pb.rand(10,10),cmap=cmap_RB)
     pb.colorbar()
diff --git a/GPy/plotting/matplot_dep/base_plots.py b/GPy/plotting/matplot_dep/base_plots.py
index 382d530d..45597628 100644
--- a/GPy/plotting/matplot_dep/base_plots.py
+++ b/GPy/plotting/matplot_dep/base_plots.py
@@ -4,7 +4,7 @@
 
 try:
     #import Tango
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 import numpy as np
diff --git a/GPy/plotting/matplot_dep/inference_plots.py b/GPy/plotting/matplot_dep/inference_plots.py
index 02007390..ab6e5eaf 100644
--- a/GPy/plotting/matplot_dep/inference_plots.py
+++ b/GPy/plotting/matplot_dep/inference_plots.py
@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 try:
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 #import numpy as np
diff --git a/GPy/plotting/matplot_dep/kernel_plots.py b/GPy/plotting/matplot_dep/kernel_plots.py
index aa015009..c8977cef 100644
--- a/GPy/plotting/matplot_dep/kernel_plots.py
+++ b/GPy/plotting/matplot_dep/kernel_plots.py
@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-import pylab as pb
+from matplotlib import pyplot as pb
 import Tango
 from matplotlib.textpath import TextPath
 from matplotlib.transforms import offset_copy
diff --git a/GPy/plotting/matplot_dep/mapping_plots.py b/GPy/plotting/matplot_dep/mapping_plots.py
index 53bc1de2..641b3fae 100644
--- a/GPy/plotting/matplot_dep/mapping_plots.py
+++ b/GPy/plotting/matplot_dep/mapping_plots.py
@@ -4,7 +4,7 @@
 import numpy as np
 try:
     import Tango
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 from base_plots import x_frame1D, x_frame2D
diff --git a/GPy/plotting/matplot_dep/maps.py b/GPy/plotting/matplot_dep/maps.py
index f38c5d93..faec9eb3 100644
--- a/GPy/plotting/matplot_dep/maps.py
+++ b/GPy/plotting/matplot_dep/maps.py
@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 import numpy as np
 try:
-    import pylab as pb
+    from matplotlib import pyplot as pb
     from matplotlib.patches import Polygon
     from matplotlib.collections import PatchCollection
     #from matplotlib import cm
diff --git a/GPy/plotting/matplot_dep/priors_plots.py b/GPy/plotting/matplot_dep/priors_plots.py
index 39dad631..51d84f86 100644
--- a/GPy/plotting/matplot_dep/priors_plots.py
+++ b/GPy/plotting/matplot_dep/priors_plots.py
@@ -4,7 +4,7 @@
 
 import numpy as np
 try:
-    import pylab as pb
+    from matplotlib import pyplot as pb
 except:
     pass
 
diff --git a/GPy/plotting/matplot_dep/svig_plots.py b/GPy/plotting/matplot_dep/svig_plots.py
index 95344643..7152e8a6 100644
--- a/GPy/plotting/matplot_dep/svig_plots.py
+++ b/GPy/plotting/matplot_dep/svig_plots.py
@@ -2,7 +2,7 @@
 # Licensed under the BSD 3-clause license (see LICENSE.txt)
 
 import numpy as np
-import pylab as pb
+from matplotlib import pyplot as pb
 
 
 def plot(model, ax=None, fignum=None, Z_height=None, **kwargs):
diff --git a/GPy/plotting/matplot_dep/variational_plots.py b/GPy/plotting/matplot_dep/variational_plots.py
index 5a62e5d8..24e613aa 100644
--- a/GPy/plotting/matplot_dep/variational_plots.py
+++ b/GPy/plotting/matplot_dep/variational_plots.py
@@ -1,4 +1,4 @@
-import pylab as pb, numpy as np
+from matplotlib import pyplot as pb, numpy as np
 
 def plot(parameterized, fignum=None, ax=None, colors=None, figsize=(12, 6)):
     """
diff --git a/GPy/testing/likelihood_tests.py b/GPy/testing/likelihood_tests.py
index 79bbdc36..1d181040 100644
--- a/GPy/testing/likelihood_tests.py
+++ b/GPy/testing/likelihood_tests.py
@@ -799,7 +799,7 @@ class LaplaceTests(unittest.TestCase):
         post_mean_approx, post_var_approx, = m2.predict(X)
 
         if debug:
-            import pylab as pb
+            from matplotlib import pyplot as pb
             pb.figure(5)
             pb.title('posterior means')
             pb.scatter(X, post_mean, c='g')
diff --git a/GPy/util/datasets.py b/GPy/util/datasets.py
index 57755ea9..f9173365 100644
--- a/GPy/util/datasets.py
+++ b/GPy/util/datasets.py
@@ -374,7 +374,7 @@ def football_data(season='1314', data_set='football_data'):
     data_resources[data_set_season]['files'] = [files]
     if not data_available(data_set_season):
         download_data(data_set_season)
-    import pylab as pb
+    from matplotlib import pyplot as pb
     for file in reversed(files):
         filename = os.path.join(data_path, data_set_season, file)
         # rewrite files removing blank rows.