diff --git a/rdt/transformers/numerical.py b/rdt/transformers/numerical.py
index d0edb29..43b522a 100644
--- a/rdt/transformers/numerical.py
+++ b/rdt/transformers/numerical.py
@@ -383,10 +383,13 @@ class ClusterBasedNormalizer(FloatFormatter):
     """Transformer for numerical data using a Bayesian Gaussian Mixture Model.
 
     This transformation takes a numerical value and transforms it using a Bayesian GMM
-    model. It generates two outputs, a discrete value which indicates the selected
+    model. It generates two outputs: a discrete value which indicates the selected
     'component' of the GMM and a continuous value which represents the normalized value
     based on the mean and std of the selected component.
 
+    The ``STD_MULTIPLIER`` attribute (default 4) is used to scale the standard deviation
+    when normalizing values, ensuring the transformed data falls within a reasonable range.
+
     Args:
         model_missing_values (bool):
             **DEPRECATED** Whether to create a new column to indicate which values were null or
@@ -463,7 +466,6 @@ class ClusterBasedNormalizer(FloatFormatter):
             n_components=self.max_clusters,
             weight_concentration_prior_type='dirichlet_process',
             weight_concentration_prior=0.001,
-            n_init=1,
             random_state=self._get_current_random_seed()
         )
 
@@ -496,8 +498,6 @@ class ClusterBasedNormalizer(FloatFormatter):
         means = self._bgm_transformer.means_.reshape((1, self.max_clusters))
 
         stds = np.sqrt(self._bgm_transformer.covariances_).reshape((1, self.max_clusters))
-        normalized_values = (data - means) / (self.STD_MULTIPLIER * stds)
-        normalized_values = normalized_values[:, self.valid_component_indicator]
         component_probs = self._bgm_transformer.predict_proba(data)
         component_probs = component_probs[:, self.valid_component_indicator]
 
@@ -510,15 +510,22 @@ class ClusterBasedNormalizer(FloatFormatter):
                 p=component_prob_t
             )
 
-        aranged = np.arange(len(data))
-        normalized = normalized_values[aranged, selected_component].reshape([-1, 1])
+        # Map valid component indices to original cluster indices
+        valid_indices = np.where(self.valid_component_indicator)[0]
+        selected_cluster = valid_indices[selected_component]
+
+        # Use broadcasting: means shape is (1, max_clusters), select per-component means
+        # Reshape to (n, 1) for proper broadcasting with data (n, 1)
+        selected_means = means[0, selected_cluster].reshape(-1, 1)  # shape: (n, 1)
+        selected_stds = stds[0, selected_cluster].reshape(-1, 1)    # shape: (n, 1)
+        normalized = (data - selected_means) / (self.STD_MULTIPLIER * selected_stds)
         normalized = np.clip(normalized, -.99, .99)
-        normalized = normalized[:, 0]
+        normalized = normalized.reshape([-1, 1])[:, 0]
         rows = [normalized, selected_component]
         if self.null_transformer and self.null_transformer.models_missing_values():
             rows.append(model_missing_values)
 
-        return np.stack(rows, axis=1)  # noqa: PD013
+        return np.stack(rows, axis=1)
 
     def _reverse_transform_helper(self, data):
         normalized = np.clip(data[:, 0], -1, 1)
@@ -546,8 +553,6 @@ class ClusterBasedNormalizer(FloatFormatter):
 
         recovered_data = self._reverse_transform_helper(data)
         if self.null_transformer and self.null_transformer.models_missing_values():
-            data = np.stack([recovered_data, data[:, -1]], axis=1)  # noqa: PD013
-        else:
-            data = recovered_data
+            data = np.stack([recovered_data, data[:, -1]], axis=1)
 
         return super()._reverse_transform(data)