diff --git a/clumpy/allocation/__pycache__/_allocator.cpython-38.pyc b/clumpy/allocation/__pycache__/_allocator.cpython-38.pyc
index 3923f9a9c32b480666aa193583d1e23ca68ddc78..534b19d9919afd5686f303c9f221671efce3b2eb 100644
Binary files a/clumpy/allocation/__pycache__/_allocator.cpython-38.pyc and b/clumpy/allocation/__pycache__/_allocator.cpython-38.pyc differ
diff --git a/clumpy/allocation/__pycache__/_unbiased.cpython-38.pyc b/clumpy/allocation/__pycache__/_unbiased.cpython-38.pyc
index 1e566c3ce682586650f2bbf770e5137057e38290..6f9ac9deccbc990d9edb82f15c86ed632713f611 100644
Binary files a/clumpy/allocation/__pycache__/_unbiased.cpython-38.pyc and b/clumpy/allocation/__pycache__/_unbiased.cpython-38.pyc differ
diff --git a/clumpy/allocation/_allocator.py b/clumpy/allocation/_allocator.py
index 28c9d4ad7a0898afd45c99be9cc969948eafed43..e06070b3f9266256b06803cc2da1b1a123225776 100644
--- a/clumpy/allocation/_allocator.py
+++ b/clumpy/allocation/_allocator.py
@@ -12,6 +12,8 @@ from ..tools._path import path_split
 from ._gart import generalized_allocation_rejection_test
 from copy import deepcopy
 
+from scipy.stats import norm
+
 class Allocator():
     """
     Allocator
@@ -67,7 +69,47 @@ class Allocator():
     #                   mask=mask)
         
     #     return(lul, proba_layer)
+    
+    def nb_monte_carlo(self,
+                       lul:LandUseLayer,
+                       tm:TransitionMatrix,
+                       features=None,
+                       mask:MaskLayer=None,
+                       alpha = 0.05,
+                       epsilon = 0.001):
+        
+        if features is None:
+            features = self.calibrator.features
+                
+        initial_state = self.calibrator.initial_state
+        final_states = self.calibrator.tpe.get_final_states()
+        
+        final_states_id = {final_state:final_states.index(final_state) for final_state in final_states}
+        P_v = np.array([tm.get(int(initial_state),
+                               int(final_state)) for final_state in final_states])
+        
+        J = lul_origin.get_J(state=initial_state,
+                      mask=mask)
+        X = lul_origin.get_X(J=J, 
+                             features=features)
         
+        X = self.calibrator.feature_selector.transform(X)
+        
+        P, final_states, P_Y = self.calibrator.tpe.transition_probabilities(
+            J=J,
+            Y=X,
+            P_v=P_v,
+            return_P_Y=True,
+            return_P_Y__v=False)
+        
+        p_alpha = norm.ppf(1-alpha/2)
+        
+        # remove initial state
+        if initial_state in final_states:
+            P = np.delete(P, list(final_states).index(initial_state), axis=1)
+        
+        return np.max(p_alpha**2/epsilon**2 * P * (1-P) / (P_Y * P.shape[0]))
+    
     def _clean_proba(self, 
                     P, 
                     final_states):
diff --git a/clumpy/allocation/_unbiased.py b/clumpy/allocation/_unbiased.py
index 61524a0d680b43e5b0f818014f2561f935e23106..8643ff225721c73393006dada3a989753fdb179c 100644
--- a/clumpy/allocation/_unbiased.py
+++ b/clumpy/allocation/_unbiased.py
@@ -42,7 +42,7 @@ class Unbiased(Allocator):
         super().__init__(calibrator=calibrator,
                          verbose=verbose,
                          verbose_heading_level=verbose_heading_level)
-
+    
     def allocate(self,
                  lul:LandUseLayer,
                  tm:TransitionMatrix,