diff --git a/pcntoolkit/model/SHASH.py b/pcntoolkit/model/SHASH.py
index 1312020e..68bc8f69 100644
--- a/pcntoolkit/model/SHASH.py
+++ b/pcntoolkit/model/SHASH.py
@@ -330,7 +330,7 @@ def dist(cls, mu, sigma, epsilon, delta, **kwargs):
 
     def logp(value, mu, sigma, epsilon, delta):
         mean = m(epsilon, delta, 1)
-        var = m(epsilon, delta, 2)
+        var = m(epsilon, delta, 2) - mean**2
         remapped_value = ((value - mu) / sigma) * np.sqrt(var) + mean
         this_S = S(remapped_value, epsilon, delta)
         this_S_sqr = np.square(this_S)
diff --git a/pcntoolkit/model/hbr.py b/pcntoolkit/model/hbr.py
index f85748dd..ca4dfc36 100644
--- a/pcntoolkit/model/hbr.py
+++ b/pcntoolkit/model/hbr.py
@@ -392,6 +392,13 @@ def estimate(self, X, y, batch_effects, **kwargs):
                     dummy_array = xarray.DataArray(data = np.zeros((len(chain), len(draw), 1)), coords = {'chain':chain, 'draw':draw,'empty':np.array([0])}, name=j)
                     self.idata.posterior[j] = dummy_array
                     self.vars_to_sample.append(j)
+
+            # zero-out all data
+            for i in self.idata.constant_data.data_vars:
+                self.idata.constant_data[i] *= 0
+            for i in self.idata.observed_data.data_vars:
+                self.idata.observed_data[i] *= 0
+
         return self.idata
 
     def predict(