diff --git a/GPy/core/model.py b/GPy/core/model.py
index 0c1594a2..46cf6ac9 100644
--- a/GPy/core/model.py
+++ b/GPy/core/model.py
@@ -308,9 +308,7 @@ class model(parameterised):
         numerical_gradient = (f1-f2)/(2*dx)
         ratio = (f1-f2)/(2*np.dot(dx,gradient))
         if verbose:
-            #print "gradient = ",gradient
-            #print "numerical gradient = ",numerical_gradient
-            print " Gradient ratio = ", ratio, '\n'
+            print "Gradient ratio = ", ratio, '\n'
             sys.stdout.flush()
 
         if (np.abs(1.-ratio)<tolerance) and not np.isnan(ratio):
@@ -319,10 +317,24 @@ class model(parameterised):
         else:
             if verbose:
                 print "Global check failed. Testing individual gradients\n"
-            try:
-                names = self.extract_param_names()
-            except NotImplementedError:
-                names = ['Variable %i'%i for i in range(len(x))]
+
+                try:
+                    names = self.extract_param_names()
+                except NotImplementedError:
+                    names = ['Variable %i'%i for i in range(len(x))]
+
+                # Prepare for pretty-printing
+                header = ['Name', 'Ratio', 'Difference', 'Analytical', 'Numerical']
+                max_names = max([len(names[i]) for i in range(len(names))] + [len(header[0])])
+                float_len = 10
+                cols = [max_names]
+                cols.extend([max(float_len, len(header[i])) for i in range(1, len(header))])
+                cols = np.array(cols) + 5
+                header_string = ["{h:^{col}}".format(h = header[i], col = cols[i]) for i in range(len(cols))]
+                header_string = map(lambda x: '|'.join(x), [header_string])
+                separator = '-'*len(header_string[0])
+                print '\n'.join([header_string[0], separator])
+
             for i in range(len(x)):
                 xx = x.copy()
                 xx[i] += step
@@ -338,13 +350,20 @@ class model(parameterised):
                 numerical_gradient = (f1-f2)/(2*step)
                 ratio = (f1-f2)/(2*step*gradient)
                 difference = np.abs((f1-f2)/2/step - gradient)
+
                 if verbose:
                     if (np.abs(ratio-1)<tolerance):
-                        formatted_name = "\033[92m {0:10s} \033[0m".format(names[i])
+                        formatted_name = "\033[92m {0} \033[0m".format(names[i])
                     else:
-                        formatted_name = "\033[91m {0:10s} \033[0m".format(names[i])
-                        
-                    print formatted_name + "   ratio: {0:15f}   difference: {1:15f}   analytical: {2:15f}    numerical: {3:15f}\n".format(float(ratio), float(difference), gradient, float(numerical_gradient)),
-        
+                        formatted_name = "\033[91m {0} \033[0m".format(names[i])
+                    r = '%.6f' % float(ratio)
+                    d = '%.6f' % float(difference)
+                    g = '%.6f' % gradient
+                    ng = '%.6f' % float(numerical_gradient)
+                    grad_string = "{0:^{c0}}|{1:^{c1}}|{2:^{c2}}|{3:^{c3}}|{4:^{c4}}".format(formatted_name,r,d,g, ng, c0 = cols[0]+9, c1 = cols[1], c2 = cols[2], c3 = cols[3], c4 = cols[4])
+                    print grad_string
+
+            print ''
+            
             return False
         return True