diff --git a/rdt/transformers/numerical.py b/rdt/transformers/numerical.py
index d0edb29..7a27367 100644
--- a/rdt/transformers/numerical.py
+++ b/rdt/transformers/numerical.py
@@ -382,10 +382,16 @@ class GaussianNormalizer(FloatFormatter):
 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
-    'component' of the GMM and a continuous value which represents the normalized value
-    based on the mean and std of the selected component.
+    During fitting, this transformer models the data using a Bayesian Gaussian Mixture
+    (BGMM). Each BGMM component represents a cluster of similar data points.
+
+    For transforming new data:
+    1. Assign each sample to a BGMM component based on posterior probabilities
+    2. Normalize the assigned component's mean and standard deviation to [-1, 1] range
+
+    The ``_transform`` method returns a 2-column array where:
+    - First column: normalized values in range [-1, 1] (float)
+    - Second column: selected component indices (int)
 
     Args:
         model_missing_values (bool):
@@ -426,6 +432,9 @@ class ClusterBasedNormalizer(FloatFormatter):
             than the ``weight_threshold``, it's indicated with True, otherwise it's set to False.
     """
 
+    # Multiplier used to convert component standard deviations to a normalized range.
+    # A value of 4 ensures that approximately 95% of values from each component
+    # fall within the [-1, 1] range after normalization.
     STD_MULTIPLIER = 4
     _bgm_transformer = None
     valid_component_indicator = None
@@ -463,7 +472,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()
         )
 
@@ -472,6 +480,7 @@ class ClusterBasedNormalizer(FloatFormatter):
         if data.ndim > 1:
             data = data[:, 0]
 
+        # Ignore warnings during BGMM fitting to avoid convergence warnings
         with warnings.catch_warnings():
             warnings.simplefilter('ignore')
             self._bgm_transformer.fit(data.reshape(-1, 1))
@@ -493,11 +502,8 @@ class ClusterBasedNormalizer(FloatFormatter):
             data, model_missing_values = data[:, 0], data[:, 1]
 
         data = data.reshape((len(data), 1))
-        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]
+        
+        # First select components based on probabilities
         component_probs = self._bgm_transformer.predict_proba(data)
         component_probs = component_probs[:, self.valid_component_indicator]
 
@@ -510,15 +516,27 @@ class ClusterBasedNormalizer(FloatFormatter):
                 p=component_prob_t
             )
 
+        # Then normalize the selected component values
+        means = self._bgm_transformer.means_.reshape((1, self.max_clusters))
+        stds = np.sqrt(self._bgm_transformer.covariances_).reshape((1, self.max_clusters))
+        
+        # Compute normalized values for all components first
+        normalized_values = (data - means) / (self.STD_MULTIPLIER * stds)
+        # Filter to only valid components
+        normalized_values = normalized_values[:, self.valid_component_indicator]
+        
+        # Get the normalized values for selected components
         aranged = np.arange(len(data))
-        normalized = normalized_values[aranged, selected_component].reshape([-1, 1])
-        normalized = np.clip(normalized, -.99, .99)
-        normalized = normalized[:, 0]
+        selected_normalized = normalized_values[aranged, selected_component]
+        
+        # Normalize the data using selected component parameters
+        normalized = np.clip(selected_normalized, -.99, .99)
+        
         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 +564,8 @@ 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
+            data = np.stack([recovered_data, data[:, -1]], axis=1)
         else:
             data = recovered_data
 
-        return super()._reverse_transform(data)
+        return super()._reverse_transform(data)
\ No newline at end of file