diff --git a/pymc3/step_methods/hmc/hmc.py b/pymc3/step_methods/hmc/hmc.py
index dd5c05ab12..6c5aa9cf55 100644
--- a/pymc3/step_methods/hmc/hmc.py
+++ b/pymc3/step_methods/hmc/hmc.py
@@ -1,5 +1,3 @@
-import math
-
 import numpy as np
 
 from ..arraystep import Competence
@@ -111,7 +109,7 @@ def _hamiltonian_step(self, start, p0, step_size):
                     'Divergence encountered, large integration error.',
                     None, state)
 
-        accept_stat = min(1, math.exp(energy_change))
+        accept_stat = min(1, np.exp(energy_change))
 
         if div_info is not None or np.random.rand() >= accept_stat:
             end = start
diff --git a/pymc3/step_methods/step_sizes.py b/pymc3/step_methods/step_sizes.py
index b063216baf..6cc3de5a87 100644
--- a/pymc3/step_methods/step_sizes.py
+++ b/pymc3/step_methods/step_sizes.py
@@ -1,5 +1,3 @@
-import math
-
 import numpy as np
 from scipy import stats
 
@@ -8,22 +6,22 @@
 
 class DualAverageAdaptation(object):
     def __init__(self, initial_step, target, gamma, k, t0):
-        self._log_step = math.log(initial_step)
+        self._log_step = np.log(initial_step)
         self._log_bar = self._log_step
         self._target = target
         self._hbar = 0.
         self._k = k
         self._t0 = t0
         self._count = 1
-        self._mu = math.log(10 * initial_step)
+        self._mu = np.log(10 * initial_step)
         self._gamma = gamma
         self._tuned_stats = []
 
     def current(self, tune):
         if tune:
-            return math.exp(self._log_step)
+            return np.exp(self._log_step)
         else:
-            return math.exp(self._log_bar)
+            return np.exp(self._log_bar)
 
     def update(self, accept_stat, tune):
         if not tune:
@@ -34,15 +32,15 @@ def update(self, accept_stat, tune):
         w = 1. / (count + t0)
         self._hbar = ((1 - w) * self._hbar + w * (self._target - accept_stat))
 
-        self._log_step = self._mu - self._hbar * math.sqrt(count) / self._gamma
+        self._log_step = self._mu - self._hbar * np.sqrt(count) / self._gamma
         mk = count ** -k
         self._log_bar = mk * self._log_step + (1 - mk) * self._log_bar
         self._count += 1
 
     def stats(self):
         return {
-            'step_size': math.exp(self._log_step),
-            'step_size_bar': math.exp(self._log_bar),
+            'step_size': np.exp(self._log_step),
+            'step_size_bar': np.exp(self._log_bar),
         }
 
     def warnings(self):