diff --git a/mofapy2/core/BayesNet.py b/mofapy2/core/BayesNet.py
index 5060b6b..63c1d76 100644
--- a/mofapy2/core/BayesNet.py
+++ b/mofapy2/core/BayesNet.py
@@ -189,7 +189,7 @@ def removeInactiveFactors(self, min_r2=None, return_idx=False):
             r2 = self.calculate_variance_explained()
 
             tmp = [
-                s.where((r2[g] > min_r2).sum(axis=0) == 0)[0]
+                np.where((r2[g] > min_r2).sum(axis=0) == 0)[0]
                 for g in range(self.dim["G"])
             ]
             drop_dic["min_r2"] = list(set.intersection(*map(set, tmp)))
@@ -197,7 +197,7 @@ def removeInactiveFactors(self, min_r2=None, return_idx=False):
                 drop_dic["min_r2"] = [np.random.choice(drop_dic["min_r2"])]
 
         # Drop the factors
-        drop = s.unique(np.concatenate(list(drop_dic.values())))
+        drop = np.unique(np.concatenate(list(drop_dic.values())))
         if len(drop) > 0:
             for node in self.nodes.keys():
                 self.nodes[node].removeFactors(drop)
@@ -401,7 +401,7 @@ def print_verbose_message(self, i):
         # print('Peak memory usage: %.2f MB' % (resource.getrusage(resource.RUSAGE_SELF).ru_maxrss / infer_platform() ))
 
         # Variance explained
-        r2 = s.asarray(self.calculate_variance_explained(total=True)).mean(axis=0)
+        r2 = np.asarray(self.calculate_variance_explained(total=True)).mean(axis=0)
         r2[r2 < 0] = 0.0
         print(
             "- Variance explained:  "
diff --git a/mofapy2/core/distributions/multivariate_gaussian.py b/mofapy2/core/distributions/multivariate_gaussian.py
index b78e299..15afc9d 100644
--- a/mofapy2/core/distributions/multivariate_gaussian.py
+++ b/mofapy2/core/distributions/multivariate_gaussian.py
@@ -141,7 +141,7 @@ def updateExpectations(self):
         E = self.params["mean"]
 
         # second moment here of the marginal components: given by E(X_n^2) = E(X_n)^2 + Var(X_n)
-        E2 = s.empty((self.dim[0], self.dim[1]))
+        E2 = np.empty((self.dim[0], self.dim[1]))
         if self.axis_cov == 1:
             for i in range(self.dim[0]):
                 E2[i, :] = E[i, :] ** 2 + np.diag(self.params["cov"][i, :, :])
@@ -348,7 +348,7 @@ def updateExpectations(self):
         # Method to calculate expectation (N,D) and variance of th marginals (N,D)
 
         # first moments
-        E = s.empty((self.dim[0], self.dim[1]))
+        E = np.empty((self.dim[0], self.dim[1]))
         if self.axis_cov == 0:
             for i in range(self.dim[1]):
                 E[:, i] = self.params["K"][i, :, :].dot(self.params["alpha"][:, i])
@@ -359,7 +359,7 @@ def updateExpectations(self):
                 )
 
         # second moment here of the marginal components: given by E(X_n^2) = E(X_n)^2 + Var(X_n)
-        E2 = s.empty((self.dim[0], self.dim[1]))
+        E2 = np.empty((self.dim[0], self.dim[1]))
         if self.axis_cov == 0:
             for i in range(self.dim[1]):
                 A = np.diag(self.params["lamb"][:, i]).dot(
diff --git a/mofapy2/core/distributions/poisson.py b/mofapy2/core/distributions/poisson.py
index f493426..004c7dd 100644
--- a/mofapy2/core/distributions/poisson.py
+++ b/mofapy2/core/distributions/poisson.py
@@ -1,6 +1,7 @@
 import numpy as np
 import scipy as s
 import scipy.stats as stats
+import math
 from .basic_distributions import Distribution
 
 from mofapy2.core.utils import *
@@ -43,7 +44,7 @@ def density(self, x):
         theta = self.params["theta"].flatten()
         x = x.flatten()
         # return np.prod (stats.poisson.pmf(x,theta) )
-        return np.prod(np.divide(theta**x * np.exp(-theta), s.misc.factorial(x)))
+        return np.prod(np.divide(theta**x * np.exp(-theta), math.factorial(x)))
 
     def loglik(self, x):
         assert x.shape == self.dim, "Problem with the dimensionalities"
@@ -51,4 +52,4 @@ def loglik(self, x):
         theta = self.params["theta"].flatten()
         x = x.flatten()
         # return np.log( np.prod (stats.poisson.pmf(x,theta) ))
-        return np.sum(x * np.log(theta) - theta - np.log(s.misc.factorial(x)))
+        return np.sum(x * np.log(theta) - theta - np.log(math.factorial(x)))
diff --git a/mofapy2/core/nodes/nongaussian_nodes.py b/mofapy2/core/nodes/nongaussian_nodes.py
index 6e3f891..6493214 100644
--- a/mofapy2/core/nodes/nongaussian_nodes.py
+++ b/mofapy2/core/nodes/nongaussian_nodes.py
@@ -197,7 +197,7 @@ def __init__(self, dim, obs, groups, params=None, E=None):
 
         # Initialise the observed data
         assert np.all(
-            s.mod(self.obs, 1) == 0
+            np.mod(self.obs, 1) == 0
         ), "Data must not contain float numbers, only integers"
         assert np.all(self.obs >= 0), "Data must not contain negative numbers"
 
diff --git a/mofapy2/simulate/simulate_mofa.py b/mofapy2/simulate/simulate_mofa.py
index dd65808..ed15578 100644
--- a/mofapy2/simulate/simulate_mofa.py
+++ b/mofapy2/simulate/simulate_mofa.py
@@ -199,7 +199,7 @@ def mask_samples(sim, perc=0.2, perc_all_views=0):
     ]
     for m in range(M):
         for g in range(G):
-            data[m][g][masked_samples[m][g], :] = s.nan
+            data[m][g][masked_samples[m][g], :] = np.nan
 
     if perc_all_views > 0:
         masked_samples.all_views = [
@@ -208,6 +208,6 @@ def mask_samples(sim, perc=0.2, perc_all_views=0):
         ]
         for m in range(len(data)):
             for g in range(G):
-                data[m][g][masked_samples.all_views[g], :] = s.nan
+                data[m][g][masked_samples.all_views[g], :] = np.nan
 
     return data