diff --git a/examples/demo_pmcmc.py b/examples/demo_pmcmc.py
index 630c24c..980cf43 100644
--- a/examples/demo_pmcmc.py
+++ b/examples/demo_pmcmc.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 import numpy as np
 import matplotlib.pyplot as plt
 from matplotlib.patches import Polygon
@@ -84,7 +85,7 @@ def sample_z_given_x(z_curr, x,
 
         # Resample transition noise
         sigmas = resample_transition_noise(prop, z_smpls[s,:,:])
-        print "Sigmas: ", sigmas
+        print("Sigmas: ", sigmas)
         prop.set_sigma(np.sqrt(sigmas))
 
 
diff --git a/hips/distributions/circular_distribution.py b/hips/distributions/circular_distribution.py
index 92b136a..8aad67a 100644
--- a/hips/distributions/circular_distribution.py
+++ b/hips/distributions/circular_distribution.py
@@ -4,6 +4,7 @@
 Scott Linderman
 2014
 """
+from __future__ import print_function
 import numpy as np
 import scipy.interpolate
 import matplotlib.patches
@@ -277,5 +278,5 @@ def test_circular_distribution():
     cmap = matplotlib.cm.get_cmap('Greys')
     cd.plot(show=True, plot_data=False, cmap=cmap)
 
-    print np.mean(cd.areas.ravel())
-    print np.std(cd.areas.ravel())
\ No newline at end of file
+    print(np.mean(cd.areas.ravel()))
+    print(np.std(cd.areas.ravel()))
\ No newline at end of file
diff --git a/hips/distributions/polya_gamma.py b/hips/distributions/polya_gamma.py
index 5336bba..d21a786 100644
--- a/hips/distributions/polya_gamma.py
+++ b/hips/distributions/polya_gamma.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 import numpy as np
 import numpy.random as npr
 
@@ -78,7 +79,7 @@ def test_polya_gamma():
     sortby = 'cumulative'
     ps = pstats.Stats(pr, stream=s).sort_stats(sortby)
     ps.print_stats()
-    print s.getvalue()
+    print(s.getvalue())
 
     # import matplotlib.pyplot as plt
     # plt.figure()
diff --git a/hips/inference/ais.py b/hips/inference/ais.py
index c8ac9ee..7cc8aa5 100644
--- a/hips/inference/ais.py
+++ b/hips/inference/ais.py
@@ -3,12 +3,15 @@
 here to be used directly since you'll need to rewrite the samplers for each model,
 but perhaps still a useful reference.
 """
+from __future__ import print_function
+from __future__ import absolute_import
 import numpy as np
 from scipy.misc import logsumexp
 import matplotlib.pyplot as plt
 
-from mh import mh
-from hmc import hmc
+from .mh import mh
+from .hmc import hmc
+from functools import reduce
 
 # Let's work with a Gaussian-Gaussian model. Say the mean of the Gaussian is
 # distributed according to a mean zero Gaussian, and the likelihood is a
@@ -70,7 +73,7 @@ def ais(eta):
 
     # Sample m points
     for m in range(M):
-        print "M: %d" % m
+        print("M: %d" % m)
         # Sample mus from each of the intermediate distributions,
         # starting with a draw from the prior.
         mus = np.zeros(N)
@@ -121,7 +124,7 @@ def ais(eta):
 for j,eta in enumerate(etas):
     # Compute true evidence. We can calculate it in closed form or by plugging in mu's.
     # Since the evidence is not a function of mu, it should be the same for all mus.
-    print 'eta: %f' % eta
+    print('eta: %f' % eta)
     mu = 0.0
     true_evidence[j] = log_prior(mu, eta) + log_lkhd(mu) - true_log_posterior(mu, eta)
 
diff --git a/hips/inference/ars.py b/hips/inference/ars.py
index 5093f28..22b099c 100644
--- a/hips/inference/ars.py
+++ b/hips/inference/ars.py
@@ -1,6 +1,7 @@
 """
 Adaptive rejection sampling
 """
+from __future__ import print_function
 import numpy as np
 from scipy.misc import logsumexp
 
@@ -94,7 +95,7 @@ def adaptive_rejection_sample(func, xs, v_xs, domain, stepsz=1.0, debug=False):
         if u <= lhVal - uhVal:
             # accept, u is below lower bound
             if debug:
-                print "Sample found after %d rejects" % rejects
+                print("Sample found after %d rejects" % rejects)
             return x
 
         # Otherwise we must compute the actual function
@@ -102,7 +103,7 @@ def adaptive_rejection_sample(func, xs, v_xs, domain, stepsz=1.0, debug=False):
         if u <= vx - uhVal:
             # accept, u is between lower bound and f
             if debug:
-                print "Sample found after %d rejects" % rejects
+                print("Sample found after %d rejects" % rejects)
             return x
 
         # If we made it this far, we rejected.
@@ -122,7 +123,7 @@ def adaptive_rejection_sample(func, xs, v_xs, domain, stepsz=1.0, debug=False):
         lowerHull, upperHull = _ars_compute_hulls(xs, v_xs, domain)
     
         if debug:
-            print 'reject %d' % rejects
+            print('reject %d' % rejects)
             
         rejects += 1
 
@@ -145,7 +146,7 @@ def _check_concavity(xs, v_xs):
     g2 = np.diff(g)/np.diff(xs[:-1])
     if not np.all(g2<=0):
         # raise Exception("It looks like the function to be sampled is not log concave!")
-        print Warning("It looks like the function to be sampled is not log concave!")
+        print(Warning("It looks like the function to be sampled is not log concave!"))
 
 def _check_boundary_grads(func, xs, v_xs, domain, stepsz):
     """
@@ -454,14 +455,14 @@ def _signed_lse(m, b, a1, a0):
     lse = b - np.log(m*sgn) + am + np.log(se*sgn)
 
     if not np.isfinite(lse):
-        print "LSE is not finite"
-        print "lse: %f" % lse
-        print "m: %f" % m
-        print "b: %f" % b
-        print "a1: %f" % a1
-        print "a2: %f" % a0
-        print "am: %f" % am
-        print "se: %f" % se
+        print("LSE is not finite")
+        print("lse: %f" % lse)
+        print("m: %f" % m)
+        print("b: %f" % b)
+        print("a1: %f" % a1)
+        print("a2: %f" % a0)
+        print("am: %f" % am)
+        print("se: %f" % se)
         raise Exception("LSE is not finite!")
     
     return lse
@@ -583,10 +584,10 @@ def _ars_compute_hulls(S, fS, domain):
         h.pr = prs[i]
 
     if not np.all(np.isfinite(prs)):
-        print "ARS prs contains Inf or NaN"
-        print lnprs
-        print lnZ
-        print prs
+        print("ARS prs contains Inf or NaN")
+        print(lnprs)
+        print(lnZ)
+        print(prs)
         import pdb; pdb.set_trace()
         raise Exception("ARS prs contains Inf or NaN")
 
@@ -595,12 +596,12 @@ def _ars_compute_hulls(S, fS, domain):
 def _ars_sample_upper_hull(upperHull):
     prs = np.array([h.pr for h in upperHull])
     if not np.all(np.isfinite(prs)):
-        print prs
+        print(prs)
         raise Exception("ARS prs contains Inf or NaN")
 
     cdf = np.cumsum(prs)
     if not np.all(np.isfinite(cdf)):
-        print cdf
+        print(cdf)
         raise Exception("ARS cumsum Inf or NaN")
     
     # randomly choose a line segment
@@ -833,7 +834,7 @@ def test_gamma_linear_regression_ars():
     y_smpls[0] = c_smpls[0]*x + sig*np.random.randn()
     for s in np.arange(1,N_samples):
         if np.mod(s, 100) == 0:
-            print "Sample ", s
+            print("Sample ", s)
         # Sample y given c
         y_smpls[s] = c_smpls[s-1]*x + sig*np.random.randn()
 
diff --git a/hips/inference/ars2.py b/hips/inference/ars2.py
index f87c73f..aa665a5 100644
--- a/hips/inference/ars2.py
+++ b/hips/inference/ars2.py
@@ -1,6 +1,7 @@
 """
 Adaptive rejection sampling
 """
+from __future__ import print_function
 import numpy as np
 from scipy.misc import logsumexp
 
@@ -91,7 +92,7 @@ def initialize_hull_points_and_domain(self, stepsz=1.0):
                         continue
                     # If the function is not finite, move our domain in
                     else:
-                        print "Found left bound at ", left_hull_pnt
+                        print("Found left bound at ", left_hull_pnt)
                         self.lb = left_hull_pnt
             left_initialized = True
 
@@ -144,7 +145,7 @@ def sample(self, debug=False):
             if u <= lhVal - uhVal:
                 # accept, u is below lower bound
                 if debug:
-                    print "Sample found after %d rejects" % rejects
+                    print("Sample found after %d rejects" % rejects)
                 return x
 
             # Otherwise we must compute the actual function
@@ -152,7 +153,7 @@ def sample(self, debug=False):
             if u <= vx - uhVal:
                 # accept, u is between lower bound and f
                 if debug:
-                    print "Sample found after %d rejects" % rejects
+                    print("Sample found after %d rejects" % rejects)
                 return x
 
             # If we made it this far, we rejected.
@@ -172,7 +173,7 @@ def sample(self, debug=False):
             self.compute_hulls()
 
             if debug:
-                print 'reject %d' % rejects
+                print('reject %d' % rejects)
 
             rejects += 1
 
@@ -330,10 +331,10 @@ def compute_hulls(self):
             h.pr = prs[i]
 
         if not np.all(np.isfinite(prs)):
-            print "ARS prs contains Inf or NaN"
-            print lnprs
-            print lnZ
-            print prs
+            print("ARS prs contains Inf or NaN")
+            print(lnprs)
+            print(lnZ)
+            print(prs)
             import pdb; pdb.set_trace()
             raise Exception("ARS prs contains Inf or NaN")
 
@@ -344,12 +345,12 @@ def sample_upper_hull(self):
         upperHull = self.upper_hull
         prs = np.array([h.pr for h in upperHull])
         if not np.all(np.isfinite(prs)):
-            print prs
+            print(prs)
             raise Exception("ARS prs contains Inf or NaN")
 
         cdf = np.cumsum(prs)
         if not np.all(np.isfinite(cdf)):
-            print cdf
+            print(cdf)
             raise Exception("ARS cumsum Inf or NaN")
 
         # randomly choose a line segment
@@ -458,14 +459,14 @@ def _signed_lse(m, b, a1, a0):
     lse = b - np.log(m*sgn) + am + np.log(se*sgn)
 
     if not np.isfinite(lse):
-        print "LSE is not finite"
-        print "lse: %f" % lse
-        print "m: %f" % m
-        print "b: %f" % b
-        print "a1: %f" % a1
-        print "a2: %f" % a0
-        print "am: %f" % am
-        print "se: %f" % se
+        print("LSE is not finite")
+        print("lse: %f" % lse)
+        print("m: %f" % m)
+        print("b: %f" % b)
+        print("a1: %f" % a1)
+        print("a2: %f" % a0)
+        print("am: %f" % am)
+        print("se: %f" % se)
         raise Exception("LSE is not finite!")
 
     return lse
diff --git a/hips/inference/elliptical_slice.py b/hips/inference/elliptical_slice.py
index c6479dc..1e30890 100644
--- a/hips/inference/elliptical_slice.py
+++ b/hips/inference/elliptical_slice.py
@@ -46,6 +46,7 @@
 Intelligence and Statistics (AISTATS), JMLR W&CP 9:541-548, 2010.
 
 """
+from __future__ import print_function
 import math
 import numpy as np
 
@@ -62,7 +63,7 @@ def elliptical_slice(xx, prior, log_like_fn, cur_log_like=None, angle_range=0, l
     else:
         # User provided Cholesky of prior covariance
         if np.shape(prior) != (D,D):
-            print "Prior shape: ", prior.shape
+            print("Prior shape: ", prior.shape)
             raise Exception("Prior must be given by a D-element sample or DxD chol(Sigma, 'lower')")
             
         nu = np.reshape(np.dot(prior, np.random.randn(D,1)).T, np.shape(xx))
diff --git a/hips/inference/hmc.py b/hips/inference/hmc.py
index c49f73e..f8640da 100644
--- a/hips/inference/hmc.py
+++ b/hips/inference/hmc.py
@@ -6,6 +6,7 @@
 slinderman@seas.harvard.edu
 2012-2014
 """
+from __future__ import print_function
 import numpy as np
 
 def hmc(U, 
@@ -149,7 +150,7 @@ def test_gamma_linear_regression_hmc():
 
     for s in np.arange(1,N_samples):
         if np.mod(s, 100) == 0:
-            print "Sample ", s
+            print("Sample ", s)
         # Sample y given c
         y_smpls[s] = np.exp(logc_smpls[s-1])*x + sig*np.random.randn()
 
diff --git a/hips/inference/slicesample.py b/hips/inference/slicesample.py
index 513205f..e744655 100644
--- a/hips/inference/slicesample.py
+++ b/hips/inference/slicesample.py
@@ -6,6 +6,7 @@
 slinderman@seas.harvard.edu
 2013-2014
 """
+from __future__ import print_function
 import numpy as np
 
 def slicesample(xx, llh_func, last_llh=None, step=1, step_out=True, x_l=None, x_r=None, lb=-np.Inf, ub=np.Inf):
@@ -88,7 +89,7 @@ def test_slicesample():
     n_iter = 1000
 
     # Gamma distribution (bounded on left)
-    print "Gamma test"
+    print("Gamma test")
     g = gamma(2.0, loc=0., scale=2.0)
 
     smpls = np.zeros(n_iter)
@@ -97,10 +98,10 @@ def test_slicesample():
         sn, _ = slicesample(smpls[n-1], g.logpdf, lb=1e-5)
         smpls[n] = sn
 
-    print "Expected gamma mean: ", g.mean()
-    print "Inferred gamma mean: ", smpls.mean()
-    print "Expected gamma std:  ", g.std()
-    print "Inferred gamma std:  ", smpls.std()
+    print("Expected gamma mean: ", g.mean())
+    print("Inferred gamma mean: ", smpls.mean())
+    print("Expected gamma std:  ", g.std())
+    print("Inferred gamma std:  ", smpls.std())
 
     fig, ax = plt.subplots(1, 1)
     x = np.linspace(1e-5, g.mean() + 4*g.std(), 1000)
diff --git a/hips/movies/examples/dynamic_spike_train.py b/hips/movies/examples/dynamic_spike_train.py
index 27429c2..b617b7e 100644
--- a/hips/movies/examples/dynamic_spike_train.py
+++ b/hips/movies/examples/dynamic_spike_train.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 import matplotlib.pyplot as plt
 import numpy as np
 from hips.movies.moviemaker import *
@@ -32,7 +33,7 @@ def dynamic_spike_train_demo():
     with writer.saving(fig, "spike_train.mp4", 200):
         ti = 1
         for fri, tfr in enumerate(tframes):
-            print "Frame %d/%d" % (fri, nframes)
+            print("Frame %d/%d" % (fri, nframes))
 
             # Update the data until we exceed the frame time
             while treallife[ti] < tfr:
diff --git a/hips/movies/moviemaker.py b/hips/movies/moviemaker.py
index 5da7b5b..c1c7962 100644
--- a/hips/movies/moviemaker.py
+++ b/hips/movies/moviemaker.py
@@ -1,3 +1,4 @@
+from __future__ import print_function
 # This example uses a MovieWriter directly to grab individual frames and
 # write them to a file. This avoids any event loop integration, but has
 # the advantage of working with even the Agg backend. This is not recommended
@@ -67,7 +68,7 @@ def update(self):
         """
         self.offset += 1
         if self.offset > self.T:
-            print "WARNING: BlinkyCircleEffect has passed the length of the data!"
+            print("WARNING: BlinkyCircleEffect has passed the length of the data!")
 
         self.p.set_array(self.data[:,self.offset])
 
@@ -123,7 +124,7 @@ def update(self):
         self.st.set_data(self._kron_matrix(self._pad_matrix(self.offset)))
 
         if self.offset > self.T:
-            print "WARNING: SlidingMatrixEffect has passed the length of the data!"
+            print("WARNING: SlidingMatrixEffect has passed the length of the data!")
 
 
 class DisappearingTraceEffect(MovieEffect):
@@ -177,7 +178,7 @@ def update(self):
             self.ls[w].set_data(trace[w,0], trace[w,1])
 
         if self.offset > self.T+self.window:
-            print "WARNING: DisappearingTraceEffect has passed the length of the data!"
+            print("WARNING: DisappearingTraceEffect has passed the length of the data!")