diff --git a/Python-module/SpatialDE/util.py b/Python-module/SpatialDE/util.py
index 3962d73..55a374f 100644
--- a/Python-module/SpatialDE/util.py
+++ b/Python-module/SpatialDE/util.py
@@ -3,61 +3,62 @@
 
 
 def qvalue(pv, pi0=None):
-    '''
-    Estimates q-values from p-values
+    """
+    Estimates q-values from p-values.
 
-    This function is modified based on https://github.com/nfusi/qvalue
+    Args:
+        pv (ndarray): Array of p-values.
+        pi0 (float): If None, it's estimated as suggested in Storey and Tibshirani, 2003.
 
-    Args
-    ====
-    pi0: if None, it's estimated as suggested in Storey and Tibshirani, 2003.
+    Returns:
+        ndarray: Array of q-values with the same shape as pv.
 
-    '''
-    assert(pv.min() >= 0 and pv.max() <= 1), "p-values should be between 0 and 1"
+    Raises:
+        ValueError: If pv is not a NumPy array or has invalid values.
+        ValueError: If pi0 is not between 0 and 1.
+
+    """
+    if not isinstance(pv, sp.ndarray):
+        raise ValueError("pv should be a NumPy array")
+
+    if not (pv.min() >= 0 and pv.max() <= 1):
+        raise ValueError("p-values should be between 0 and 1")
 
     original_shape = pv.shape
-    pv = pv.ravel()  # flattens the array in place, more efficient than flatten()
+    pv = pv.ravel()
 
     m = float(len(pv))
 
-    # if the number of hypotheses is small, just set pi0 to 1
     if len(pv) < 100 and pi0 is None:
         pi0 = 1.0
     elif pi0 is not None:
         pi0 = pi0
     else:
-        # evaluate pi0 for different lambdas
-        pi0 = []
-        lam = sp.arange(0, 0.90, 0.01)
         counts = sp.array([(pv > i).sum() for i in sp.arange(0, 0.9, 0.01)])
-        for l in range(len(lam)):
-            pi0.append(counts[l]/(m*(1-lam[l])))
-
-        pi0 = sp.array(pi0)
+        lam = sp.arange(0, 0.90, 0.01)
+        pi0 = counts / (m * (1 - lam))
 
-        # fit natural cubic spline
         tck = interpolate.splrep(lam, pi0, k=3)
         pi0 = interpolate.splev(lam[-1], tck)
 
         if pi0 > 1:
             pi0 = 1.0
 
-    assert(pi0 >= 0 and pi0 <= 1), "pi0 is not between 0 and 1: %f" % pi0
+    if not (0 <= pi0 <= 1):
+        raise ValueError("pi0 is not between 0 and 1: %f" % pi0)
 
     p_ordered = sp.argsort(pv)
     pv = pv[p_ordered]
-    qv = pi0 * m/len(pv) * pv
+    qv = pi0 * m / len(pv) * pv
     qv[-1] = min(qv[-1], 1.0)
 
-    for i in range(len(pv)-2, -1, -1):
-        qv[i] = min(pi0*m*pv[i]/(i+1.0), qv[i+1])
+    i = sp.arange(len(pv) - 2, -1, -1)
+    qv[i] = sp.minimum(pi0 * m * pv[i] / (i + 1.0), qv[i + 1])
 
-    # reorder qvalues
     qv_temp = qv.copy()
     qv = sp.zeros_like(qv)
     qv[p_ordered] = qv_temp
 
-    # reshape qvalues
     qv = qv.reshape(original_shape)
 
     return qv