use systematic sampling

Author: aloctavodia

pymc/smc/smc.py

@@ -274,9 +274,7 @@ def update_beta_and_weights(self):
 
     def resample(self):
         """Resample particles based on importance weights"""
-        self.resampling_indexes = self.rng.choice(
-            np.arange(self.draws), size=self.draws, p=self.weights
-        )
+        self.resampling_indexes = systematic(self.weights)
 
         self.tempered_posterior = self.tempered_posterior[self.resampling_indexes]
         self.prior_logp = self.prior_logp[self.resampling_indexes]
@@ -546,6 +544,36 @@ def sample_settings(self):
         return stats
 
 
+def systematic(weights):
+    """
+    Systematic resampling.
+
+    Parameters
+    ----------
+    weights :
+        The weights should be probabilities and the total sum should be 1.
+
+    Returns
+    -------
+    new_indices: array
+        A vector of indices in the interval 0, ..., len(normalized_weights)
+    """
+    lnw = len(weights)
+    arange = np.arange(lnw)
+    uniform = (np.random.rand(1) + arange) / lnw
+
+    idx = 0
+    weight_accu = weights[0]
+    new_indices = np.empty(lnw, dtype=int)
+    for i in arange:
+        while uniform[i] > weight_accu:
+            idx += 1
+            weight_accu += weights[idx]
+        new_indices[i] = idx
+
+    return new_indices
+
+
 def _logp_forw(point, out_vars, in_vars, shared):
     """Compile Aesara function of the model and the input and output variables.