diff --git a/python/cuml/cuml/ensemble/_gpu_tree.py b/python/cuml/cuml/ensemble/_gpu_tree.py
new file mode 100644
index 0000000000..c1132dc2ae
--- /dev/null
+++ b/python/cuml/cuml/ensemble/_gpu_tree.py
@@ -0,0 +1,418 @@
+# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION.
+# SPDX-License-Identifier: Apache-2.0
+
+"""GPU-backed proxy objects that mimic sklearn DecisionTree estimators.
+
+These are exposed via ``RandomForest{Classifier,Regressor}.estimators_``
+and route inference through FIL while exposing the full tree structure
+(topology, node values, impurity) extracted from the treelite model.
+"""
+
+import cupy as cp
+import numpy as np
+import treelite
+from scipy.sparse import csr_matrix
+
+from cuml.common.classification import decode_labels
+from cuml.internals.base import Base
+from cuml.internals.mixins import ClassifierMixin, RegressorMixin
+from cuml.internals.outputs import (
+    exit_internal_context,
+    reflect,
+    run_in_internal_context,
+)
+
+
+def _compute_max_depth(children_left, children_right):
+    """Compute the max depth of a tree from its children arrays."""
+    n_nodes = len(children_left)
+    if n_nodes == 0:
+        return 0
+    depths = np.zeros(n_nodes, dtype=np.intp)
+    for i in range(n_nodes):
+        left = children_left[i]
+        if left != -1:
+            depths[left] = depths[i] + 1
+            depths[children_right[i]] = depths[i] + 1
+    return int(depths.max())
+
+
+def _extract_tree_from_treelite(tl_model, tree_idx, n_classes, is_classifier):
+    """Extract sklearn-compatible tree arrays from a treelite tree accessor.
+
+    Returns a dict with keys matching sklearn Tree attributes.
+    """
+    ta = tl_model.get_tree_accessor(tree_idx)
+    header = tl_model.get_header_accessor()
+
+    node_type = ta.get_field("node_type")  # 1=internal, 0=leaf
+    n_nodes = int(ta.get_field("num_nodes")[0])
+    cleft = ta.get_field("cleft").astype(np.intp)
+    cright = ta.get_field("cright").astype(np.intp)
+    split_index = ta.get_field("split_index").astype(np.intp)
+    threshold_raw = ta.get_field("threshold")
+    data_count = ta.get_field("data_count").astype(np.intp)
+    leaf_vector = ta.get_field("leaf_vector")
+
+    n_features = int(header.get_field("num_feature")[0])
+    is_leaf = node_type == 0
+
+    # Convert to sklearn conventions
+    children_left = cleft.copy()
+    children_right = cright.copy()
+    children_left[is_leaf] = -1
+    children_right[is_leaf] = -1
+
+    feature = split_index.copy()
+    feature[is_leaf] = -2
+
+    threshold = threshold_raw.copy().astype(np.float64)
+    threshold[is_leaf] = -2.0
+
+    n_node_samples = data_count.copy()
+
+    n_outputs = 1
+    if is_classifier:
+        leaf_probs = leaf_vector.reshape(-1, n_classes)
+
+        value = np.zeros((n_nodes, n_outputs, n_classes), dtype=np.float64)
+        leaf_counter = 0
+        for i in range(n_nodes):
+            if is_leaf[i]:
+                value[i, 0, :] = leaf_probs[leaf_counter] * data_count[i]
+                leaf_counter += 1
+
+        for i in range(n_nodes - 1, -1, -1):
+            if not is_leaf[i]:
+                value[i, 0, :] = value[cleft[i], 0, :] + value[cright[i], 0, :]
+    else:
+        n_classes = 1
+        leaf_vals = ta.get_field("leaf_value").astype(np.float64)
+
+        value = np.zeros((n_nodes, n_outputs, 1), dtype=np.float64)
+        for i in range(n_nodes):
+            if is_leaf[i]:
+                value[i, 0, 0] = leaf_vals[i] * data_count[i]
+
+        for i in range(n_nodes - 1, -1, -1):
+            if not is_leaf[i]:
+                value[i, 0, 0] = value[cleft[i], 0, 0] + value[cright[i], 0, 0]
+
+    # Impurity: Gini for classification, variance for regression
+    impurity = np.zeros(n_nodes, dtype=np.float64)
+    if is_classifier:
+        for i in range(n_nodes):
+            total = value[i, 0, :].sum()
+            if total > 0:
+                p = value[i, 0, :] / total
+                impurity[i] = 1.0 - np.sum(p**2)
+    else:
+        pass  # TODO: compute variance-based impurity for regression
+
+    max_depth = _compute_max_depth(children_left, children_right)
+
+    return {
+        "children_left": children_left,
+        "children_right": children_right,
+        "feature": feature,
+        "threshold": threshold,
+        "value": value,
+        "impurity": impurity,
+        "n_node_samples": n_node_samples,
+        "weighted_n_node_samples": n_node_samples.astype(np.float64),
+        "node_count": n_nodes,
+        "n_features": n_features,
+        "n_classes": np.array(
+            [n_classes] if is_classifier else [1], dtype=np.intp
+        ),
+        "n_outputs": n_outputs,
+        "max_depth": max_depth,
+    }
+
+
+class GPUTree:
+    """Duck-typed proxy for ``sklearn.tree._tree.Tree``.
+
+    Exposes tree structure arrays and routes ``predict``/``apply`` through FIL.
+    """
+
+    def __init__(self, tree_data, fil_model, tree_idx, parent_estimator):
+        self._fil_model = fil_model
+        self._tree_idx = tree_idx
+        self._parent_estimator = parent_estimator
+
+        self.children_left = tree_data["children_left"]
+        self.children_right = tree_data["children_right"]
+        self.feature = tree_data["feature"]
+        self.threshold = tree_data["threshold"]
+        self.value = tree_data["value"]
+        self.impurity = tree_data["impurity"]
+        self.n_node_samples = tree_data["n_node_samples"]
+        self.weighted_n_node_samples = tree_data["weighted_n_node_samples"]
+        self.node_count = tree_data["node_count"]
+        self.capacity = tree_data["node_count"]
+        self.n_features = tree_data["n_features"]
+        self.n_classes = tree_data["n_classes"]
+        self.n_outputs = tree_data["n_outputs"]
+        self.max_depth = tree_data["max_depth"]
+        self.max_n_classes = int(self.n_classes.max())
+
+    @property
+    def n_leaves(self):
+        return int(np.sum(self.children_left == -1))
+
+    def _coerce_output(self, gpu_array, X):
+        """Convert a cupy array to the appropriate output type based on the parent estimator."""
+        output_type = self._parent_estimator._get_output_type(X)
+        if output_type in (None, "cupy"):
+            return gpu_array
+        return cp.asnumpy(gpu_array)
+
+    def predict(self, X):
+        """Return per-class probabilities via FIL predict_per_tree."""
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        per_tree = self._fil_model.predict_per_tree(X_gpu)
+        result = per_tree[:, self._tree_idx]
+        if result.ndim == 1:
+            result = result[:, np.newaxis]
+        return self._coerce_output(result, X)
+
+    def apply(self, X):
+        """Return leaf node IDs via FIL apply."""
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        leaf_ids = self._fil_model.apply(X_gpu)
+        result = leaf_ids[:, self._tree_idx].astype(cp.intp)
+        return self._coerce_output(result, X)
+
+    def decision_path(self, X):
+        """CPU fallback: traverse the stored tree structure."""
+        X = np.asarray(X, dtype=np.float32)
+        n_samples = X.shape[0]
+        indptr = [0]
+        indices = []
+
+        for sample_idx in range(n_samples):
+            node = 0
+            while node != -1:
+                indices.append(node)
+                if self.children_left[node] == -1:
+                    break
+                if X[sample_idx, self.feature[node]] <= self.threshold[node]:
+                    node = self.children_left[node]
+                else:
+                    node = self.children_right[node]
+            indptr.append(len(indices))
+
+        return csr_matrix(
+            (np.ones(len(indices), dtype=np.uint8), indices, indptr),
+            shape=(n_samples, self.node_count),
+        )
+
+    def compute_feature_importances(self, normalize=True):
+        importances = np.zeros(self.n_features, dtype=np.float64)
+        is_leaf = self.children_left == -1
+        for i in range(self.node_count):
+            if is_leaf[i]:
+                continue
+            left = self.children_left[i]
+            right = self.children_right[i]
+            w = self.weighted_n_node_samples
+            importances[self.feature[i]] += (
+                w[i] * self.impurity[i]
+                - w[left] * self.impurity[left]
+                - w[right] * self.impurity[right]
+            )
+        if normalize:
+            total = importances.sum()
+            if total > 0:
+                importances /= total
+        return importances
+
+
+def _build_gpu_estimators(rf_model):
+    """Build list of GPU-backed DecisionTree proxy objects from a fitted RF."""
+    is_classifier = rf_model._estimator_type == "classifier"
+    tl_model = treelite.Model.deserialize_bytes(rf_model._treelite_model_bytes)
+    fil_model = rf_model._get_inference_nvforest_model()
+
+    header = tl_model.get_header_accessor()
+    n_classes = int(header.get_field("num_class")[0])
+    n_features = int(header.get_field("num_feature")[0])
+    n_trees = tl_model.num_tree
+
+    common_params = dict(
+        split_criterion=rf_model.split_criterion,
+        max_depth=rf_model.max_depth,
+        min_samples_split=rf_model.min_samples_split,
+        min_samples_leaf=rf_model.min_samples_leaf,
+        max_features=rf_model.max_features,
+        max_leaves=rf_model.max_leaves,
+        min_impurity_decrease=rf_model.min_impurity_decrease,
+        random_state=rf_model.random_state,
+    )
+
+    estimators = []
+    for tree_idx in range(n_trees):
+        tree_data = _extract_tree_from_treelite(
+            tl_model,
+            tree_idx,
+            n_classes,
+            is_classifier,
+        )
+        if is_classifier:
+            est = GPUDecisionTreeClassifier(**common_params)
+            est.classes_ = rf_model.classes_
+            est.n_classes_ = rf_model.n_classes_
+        else:
+            est = GPUDecisionTreeRegressor(**common_params)
+
+        est.tree_ = GPUTree(
+            tree_data, fil_model, tree_idx, parent_estimator=est
+        )
+        est.n_features_in_ = n_features
+        est.n_outputs_ = 1
+        est.max_features_ = n_features
+        est._fil_model = fil_model
+        est._tree_idx = tree_idx
+        estimators.append(est)
+
+    return estimators
+
+
+class GPUDecisionTreeClassifier(Base, ClassifierMixin):
+    """GPU-backed proxy for a single decision tree classifier in a forest.
+
+    This class is used to represent a fitted decision tree in the estimators_
+    attribute of a RandomForestClassifier. The goal is to make the internal
+    structure of the forest accessible. It is not a fully functional estimator.
+    For example, it does not have a `fit` method.
+    """
+
+    def __init__(
+        self,
+        split_criterion="gini",
+        max_depth=None,
+        min_samples_split=2,
+        min_samples_leaf=1,
+        max_features=None,
+        max_leaves=-1,
+        min_impurity_decrease=0.0,
+        random_state=None,
+        verbose=False,
+        output_type=None,
+    ):
+        super().__init__(verbose=verbose, output_type=output_type)
+        self.split_criterion = split_criterion
+        self.max_depth = max_depth
+        self.min_samples_split = min_samples_split
+        self.min_samples_leaf = min_samples_leaf
+        self.max_features = max_features
+        self.max_leaves = max_leaves
+        self.min_impurity_decrease = min_impurity_decrease
+        self.random_state = random_state
+
+    @classmethod
+    def _get_param_names(cls):
+        return [
+            *super()._get_param_names(),
+            "split_criterion",
+            "max_depth",
+            "min_samples_split",
+            "min_samples_leaf",
+            "max_features",
+            "max_leaves",
+            "min_impurity_decrease",
+            "random_state",
+        ]
+
+    @run_in_internal_context
+    def predict(self, X):
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        per_tree = self._fil_model.predict_per_tree(X_gpu)
+        proba = per_tree[:, self._tree_idx]
+        inds = cp.argmax(proba, axis=1)
+        with exit_internal_context():
+            output_type = self._get_output_type(X)
+        return decode_labels(inds, self.classes_, output_type=output_type)
+
+    @reflect
+    def predict_proba(self, X):
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        per_tree = self._fil_model.predict_per_tree(X_gpu)
+        return per_tree[:, self._tree_idx]
+
+    @reflect
+    def apply(self, X):
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        leaf_ids = self._fil_model.apply(X_gpu)
+        return leaf_ids[:, self._tree_idx]
+
+    @property
+    def feature_importances_(self):
+        return self.tree_.compute_feature_importances()
+
+
+class GPUDecisionTreeRegressor(Base, RegressorMixin):
+    """GPU-backed proxy for a single decision tree regressor in a forest.
+
+    This class is used to represent a fitted decision tree in the estimators_
+    attribute of a RandomForestRegressor. The goal is to make the internal
+    structure of the forest accessible. It is not a fully functional estimator.
+    For example, it does not have a `fit` method.
+    """
+
+    def __init__(
+        self,
+        split_criterion="mse",
+        max_depth=None,
+        min_samples_split=2,
+        min_samples_leaf=1,
+        max_features=None,
+        max_leaves=-1,
+        min_impurity_decrease=0.0,
+        random_state=None,
+        verbose=False,
+        output_type=None,
+    ):
+        super().__init__(verbose=verbose, output_type=output_type)
+        self.split_criterion = split_criterion
+        self.max_depth = max_depth
+        self.min_samples_split = min_samples_split
+        self.min_samples_leaf = min_samples_leaf
+        self.max_features = max_features
+        self.max_leaves = max_leaves
+        self.min_impurity_decrease = min_impurity_decrease
+        self.random_state = random_state
+
+    @classmethod
+    def _get_param_names(cls):
+        return [
+            *super()._get_param_names(),
+            "split_criterion",
+            "max_depth",
+            "min_samples_split",
+            "min_samples_leaf",
+            "max_features",
+            "max_leaves",
+            "min_impurity_decrease",
+            "random_state",
+        ]
+
+    @reflect
+    def predict(self, X):
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        per_tree = self._fil_model.predict_per_tree(X_gpu)
+        result = per_tree[:, self._tree_idx]
+        if result.ndim > 1:
+            result = result.squeeze(axis=-1)
+        return result
+
+    @reflect
+    def apply(self, X):
+        X_gpu = cp.asarray(X, dtype=cp.float32)
+        leaf_ids = self._fil_model.apply(X_gpu)
+        return leaf_ids[:, self._tree_idx]
+
+    @property
+    def feature_importances_(self):
+        return self.tree_.compute_feature_importances()
diff --git a/python/cuml/cuml/ensemble/randomforest_common.pyx b/python/cuml/cuml/ensemble/randomforest_common.pyx
index 99f5e9142d..451f6fc8df 100644
--- a/python/cuml/cuml/ensemble/randomforest_common.pyx
+++ b/python/cuml/cuml/ensemble/randomforest_common.pyx
@@ -26,6 +26,7 @@ from cuml.metrics import accuracy_score, r2_score
 from libc.stdint cimport uint64_t, uintptr_t
 from libcpp cimport bool
 from pylibraft.common.handle cimport handle_t
+
 import nvforest
 
 from cuml.internals.logger cimport level_enum
@@ -35,6 +36,8 @@ from cuml.internals.treelite cimport (
     TreeliteSerializeModelToBytes,
 )
 
+from cuml.ensemble._gpu_tree import _build_gpu_estimators
+
 
 cdef extern from "cuml/ensemble/randomforest.hpp" namespace "ML" nogil:
     cdef enum CRITERION:
@@ -345,8 +348,30 @@ class BaseRandomForestModel(Base, InteropMixin):
         self.n_streams = n_streams
         self.oob_score = oob_score
 
+    @property
+    def estimators_(self):
+        """List of GPU-backed DecisionTree proxy objects for each tree.
+
+        Each estimator exposes the same API as sklearn's
+        ``DecisionTreeClassifier`` / ``DecisionTreeRegressor`` and routes
+        inference through FIL on the GPU.
+
+        The list is constructed lazily on first access and cached.
+        """
+        if not hasattr(self, "_treelite_model_bytes"):
+            raise AttributeError(
+                f"'{type(self).__name__}' object has no attribute 'estimators_'"
+            )
+        if (cached := getattr(self, "_estimators_cache", None)) is not None:
+            return cached
+        self._estimators_cache = _build_gpu_estimators(self)
+        return self._estimators_cache
+
     def __getstate__(self):
         state = self.__dict__.copy()
+        # FIL model isn't currently pickleable
+        state.pop("_fil_model", None)
+        state.pop("_estimators_cache", None)
         # nvForest model isn't currently pickleable
         state.pop("_nvforest_model", None)
         return state
@@ -643,6 +668,9 @@ class BaseRandomForestModel(Base, InteropMixin):
         )
         self.n_outputs_ = 1
         self._treelite_model_bytes = <bytes>(tl_bytes[:tl_bytes_len])
+        # Ensure cached models are reset
+        self._fil_model = None
+        self._estimators_cache = None
         # Reload nvforest model
         self._nvforest_model = self.as_nvforest()
 
diff --git a/python/cuml/tests/test_rf_estimators.py b/python/cuml/tests/test_rf_estimators.py
new file mode 100644
index 0000000000..54ca8ecb4f
--- /dev/null
+++ b/python/cuml/tests/test_rf_estimators.py
@@ -0,0 +1,520 @@
+# SPDX-FileCopyrightText: Copyright (c) 2026, NVIDIA CORPORATION.
+# SPDX-License-Identifier: Apache-2.0
+
+"""Tests for RandomForest.estimators_ GPU-backed proxy objects."""
+
+import cupy as cp
+import numpy as np
+import pytest
+from sklearn.datasets import make_classification, make_regression
+from sklearn.ensemble import RandomForestClassifier as skRFC
+from sklearn.ensemble import RandomForestRegressor as skRFR
+
+from cuml.ensemble import RandomForestClassifier as cuRFC
+from cuml.ensemble import RandomForestRegressor as cuRFR
+from cuml.internals.base import Base
+
+pytestmark = pytest.mark.filterwarnings(
+    "ignore:The default value of 'max_depth':FutureWarning"
+)
+
+
+@pytest.fixture(scope="module")
+def clf_data():
+    X, y = make_classification(
+        n_samples=500,
+        n_features=10,
+        n_informative=5,
+        n_classes=2,
+        random_state=42,
+    )
+    return X.astype(np.float32), y.astype(np.int32)
+
+
+@pytest.fixture(scope="module")
+def multiclass_data():
+    X, y = make_classification(
+        n_samples=500,
+        n_features=10,
+        n_informative=5,
+        n_classes=5,
+        n_clusters_per_class=1,
+        random_state=42,
+    )
+    return X.astype(np.float32), y.astype(np.int32)
+
+
+@pytest.fixture(scope="module")
+def reg_data():
+    X, y = make_regression(
+        n_samples=500,
+        n_features=10,
+        n_informative=5,
+        random_state=42,
+    )
+    return X.astype(np.float32), y.astype(np.float32)
+
+
+# ---------------------------------------------------------------------------
+# GPUTree structural attributes
+# ---------------------------------------------------------------------------
+
+
+def test_tree_topology_shapes(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+
+    assert tree.node_count > 0
+    assert tree.n_features == X.shape[1]
+    assert tree.n_outputs == 1
+    assert tree.children_left.shape == (tree.node_count,)
+    assert tree.children_right.shape == (tree.node_count,)
+    assert tree.feature.shape == (tree.node_count,)
+    assert tree.threshold.shape == (tree.node_count,)
+
+
+def test_leaf_internal_consistency(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    is_leaf = tree.children_left == -1
+
+    np.testing.assert_array_equal(
+        tree.children_left[is_leaf], tree.children_right[is_leaf]
+    )
+    assert np.all(tree.feature[is_leaf] == -2)
+    assert np.all(tree.feature[~is_leaf] >= 0)
+    assert np.all(tree.children_left[~is_leaf] >= 0)
+    assert np.all(tree.children_right[~is_leaf] >= 0)
+
+
+def test_n_node_samples_parent_equals_children(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    is_leaf = tree.children_left == -1
+
+    assert tree.n_node_samples[0] > 0
+    for i in range(tree.node_count):
+        if not is_leaf[i]:
+            left = tree.children_left[i]
+            right = tree.children_right[i]
+            assert tree.n_node_samples[i] == (
+                tree.n_node_samples[left] + tree.n_node_samples[right]
+            )
+
+
+def test_estimators_count(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=7, max_depth=4, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+    assert len(rf.estimators_) == 7
+
+
+def test_estimators_attribute_error_before_fit():
+    rf = cuRFC(n_estimators=3, max_depth=4, random_state=42)
+    with pytest.raises(AttributeError, match="no attribute"):
+        rf.estimators_
+
+
+def test_max_depth_respected(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=4, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+    for est in rf.estimators_:
+        assert est.tree_.max_depth <= 4
+
+
+# ---------------------------------------------------------------------------
+# GPUTree predict / apply
+# ---------------------------------------------------------------------------
+
+
+def test_tree_predict_matches_predict_proba(clf_data):
+    """tree_.predict() returns same probabilities as estimator predict_proba."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+    tree_pred = est.tree_.predict(X)
+    est_proba = est.predict_proba(X)
+    np.testing.assert_allclose(tree_pred, est_proba, rtol=1e-5)
+
+
+def test_tree_apply_returns_valid_leaf_ids(clf_data):
+    """tree_.apply() returns leaf node IDs that are actually leaves."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    leaf_ids = tree.apply(X)
+    is_leaf = tree.children_left == -1
+
+    assert leaf_ids.shape == (X.shape[0],)
+    assert np.all(is_leaf[leaf_ids])
+
+
+def test_tree_apply_matches_estimator_apply(clf_data):
+    """tree_.apply() matches estimator-level apply()."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+    np.testing.assert_array_equal(est.tree_.apply(X), est.apply(X))
+
+
+def test_tree_predict_multiclass(multiclass_data):
+    """tree_.predict() works for multi-class problems."""
+    X, y = multiclass_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree_pred = rf.estimators_[0].tree_.predict(X[:10])
+    assert tree_pred.shape == (10, 5)
+    np.testing.assert_allclose(tree_pred.sum(axis=1), 1.0, rtol=1e-5)
+
+
+# ---------------------------------------------------------------------------
+# Node value reconstruction and impurity
+# ---------------------------------------------------------------------------
+
+
+def test_value_parent_equals_children_sum(clf_data):
+    """value[parent] == value[left] + value[right] for all internal nodes."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    for i in range(tree.node_count):
+        if tree.children_left[i] != -1:
+            left = tree.children_left[i]
+            right = tree.children_right[i]
+            np.testing.assert_allclose(
+                tree.value[i], tree.value[left] + tree.value[right], rtol=1e-5
+            )
+
+
+def test_value_sum_equals_n_node_samples(clf_data):
+    """For classification: sum(value[i]) == n_node_samples[i]."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    for est in rf.estimators_:
+        tree = est.tree_
+        for i in range(tree.node_count):
+            np.testing.assert_allclose(
+                tree.value[i].sum(),
+                tree.n_node_samples[i],
+                rtol=1e-5,
+                err_msg=f"node {i}",
+            )
+
+
+def test_impurity_bounds(clf_data):
+    """Impurity in [0, 1] for all nodes; pure leaves have impurity 0."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    assert np.all(tree.impurity >= 0)
+    assert np.all(tree.impurity <= 1)
+
+    # Pure leaves (single class) should have impurity == 0
+    is_leaf = tree.children_left == -1
+    for i in range(tree.node_count):
+        if is_leaf[i]:
+            v = tree.value[i, 0, :]
+            if np.count_nonzero(v) <= 1:
+                assert tree.impurity[i] == 0.0
+
+
+def test_from_sklearn_roundtrip_values(clf_data):
+    """from_sklearn roundtrip: cuml estimators_ predictions match sklearn's."""
+    X, y = clf_data
+    sk_rf = skRFC(n_estimators=5, max_depth=5, random_state=42)
+    sk_rf.fit(X, y)
+
+    cu_rf = cuRFC.from_sklearn(sk_rf)
+    for i in range(5):
+        sk_proba = sk_rf.estimators_[i].predict_proba(X[:20])
+        cu_proba = cu_rf.estimators_[i].predict_proba(X[:20])
+        np.testing.assert_allclose(cu_proba, sk_proba, rtol=1e-4)
+
+
+# ---------------------------------------------------------------------------
+# GPUDecisionTreeClassifier full API
+# ---------------------------------------------------------------------------
+
+
+def test_isinstance_classifier(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+
+    assert isinstance(est, Base)
+    assert hasattr(est, "tree_")
+    assert hasattr(est, "classes_")
+    assert hasattr(est, "n_classes_")
+    assert hasattr(est, "n_features_in_")
+    assert hasattr(est, "n_outputs_")
+
+
+def test_estimator_repr(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    r = repr(rf.estimators_[0])
+    assert "GPUDecisionTreeClassifier" in r
+    assert "split_criterion" in r
+
+
+def test_classifier_predict_returns_classes(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+    pred = est.predict(X)
+    assert pred.shape == (X.shape[0],)
+    assert set(pred).issubset(set(cp.asnumpy(rf.classes_)))
+
+
+def test_classifier_predict_proba_shape_and_sum(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    proba = rf.estimators_[0].predict_proba(X)
+    assert proba.shape == (X.shape[0], 2)
+    np.testing.assert_allclose(proba.sum(axis=1), 1.0, rtol=1e-5)
+
+
+def test_classifier_self_consistency(clf_data):
+    """Forest predict_proba == mean of individual estimator predict_probas."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=10, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    forest_proba = cp.asnumpy(rf.predict_proba(cp.asarray(X)))
+    est_probas = np.array([est.predict_proba(X) for est in rf.estimators_])
+    mean_proba = est_probas.mean(axis=0)
+    np.testing.assert_allclose(forest_proba, mean_proba, rtol=1e-4)
+
+
+def test_classifier_feature_importances(clf_data):
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    imp = rf.estimators_[0].feature_importances_
+    assert imp.shape == (X.shape[1],)
+    assert np.all(imp >= 0)
+    np.testing.assert_allclose(imp.sum(), 1.0, rtol=1e-5)
+
+
+def test_multiclass_predict_proba(multiclass_data):
+    X, y = multiclass_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    proba = rf.estimators_[0].predict_proba(X)
+    assert proba.shape == (X.shape[0], 5)
+    np.testing.assert_allclose(proba.sum(axis=1), 1.0, rtol=1e-5)
+
+
+# ---------------------------------------------------------------------------
+# GPUDecisionTreeRegressor
+# ---------------------------------------------------------------------------
+
+
+def test_isinstance_regressor(reg_data):
+    X, y = reg_data
+    rf = cuRFR(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+
+    assert isinstance(est, Base)
+    assert hasattr(est, "tree_")
+    assert hasattr(est, "n_features_in_")
+
+
+def test_regressor_predict_shape(reg_data):
+    X, y = reg_data
+    rf = cuRFR(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    pred = rf.estimators_[0].predict(X)
+    assert pred.shape == (X.shape[0],)
+    assert pred.dtype in (np.float32, np.float64)
+
+
+def test_regressor_self_consistency(reg_data):
+    """Forest predict == mean of individual estimator predictions."""
+    X, y = reg_data
+    rf = cuRFR(n_estimators=10, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    forest_pred = cp.asnumpy(rf.predict(cp.asarray(X))).ravel()
+    est_preds = np.array([est.predict(X) for est in rf.estimators_])
+    mean_pred = est_preds.mean(axis=0)
+    np.testing.assert_allclose(forest_pred, mean_pred, rtol=1e-4)
+
+
+def test_regressor_from_sklearn_roundtrip(reg_data):
+    X, y = reg_data
+    sk_rf = skRFR(n_estimators=5, max_depth=5, random_state=42)
+    sk_rf.fit(X, y)
+
+    cu_rf = cuRFR.from_sklearn(sk_rf)
+    for i in range(5):
+        sk_pred = sk_rf.estimators_[i].predict(X[:20])
+        cu_pred = cu_rf.estimators_[i].predict(X[:20])
+        np.testing.assert_allclose(cu_pred, sk_pred, rtol=1e-4)
+
+
+def test_regressor_tree_structure(reg_data):
+    X, y = reg_data
+    rf = cuRFR(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    is_leaf = tree.children_left == -1
+
+    assert tree.node_count > 0
+    assert np.all(tree.feature[is_leaf] == -2)
+    assert np.all(tree.feature[~is_leaf] >= 0)
+    assert tree.n_node_samples[0] > 0
+
+
+# ---------------------------------------------------------------------------
+# sklearn source-of-truth tests
+# ---------------------------------------------------------------------------
+
+
+def test_as_sklearn_roundtrip_classifier(clf_data):
+    """as_sklearn() roundtrip: exported estimators_ match our proxy."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=5, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    sk_model = rf.as_sklearn()
+    for i in range(5):
+        sk_proba = sk_model.estimators_[i].predict_proba(X[:20])
+        cu_proba = rf.estimators_[i].predict_proba(X[:20])
+        np.testing.assert_allclose(cu_proba, sk_proba, rtol=1e-4)
+
+
+def test_as_sklearn_roundtrip_regressor(reg_data):
+    """as_sklearn() roundtrip: exported estimators_ match our proxy."""
+    X, y = reg_data
+    rf = cuRFR(n_estimators=5, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    sk_model = rf.as_sklearn()
+    for i in range(5):
+        sk_pred = sk_model.estimators_[i].predict(X[:20])
+        cu_pred = rf.estimators_[i].predict(X[:20])
+        np.testing.assert_allclose(cu_pred, sk_pred, rtol=1e-4)
+
+
+def test_from_sklearn_roundtrip_classifier_apply(clf_data):
+    """from_sklearn roundtrip: apply() matches sklearn's tree."""
+    X, y = clf_data
+    sk_rf = skRFC(n_estimators=5, max_depth=5, random_state=42)
+    sk_rf.fit(X, y)
+
+    cu_rf = cuRFC.from_sklearn(sk_rf)
+    for i in range(5):
+        sk_leaf = sk_rf.estimators_[i].apply(X[:20])
+        cu_leaf = cu_rf.estimators_[i].apply(X[:20])
+        np.testing.assert_array_equal(cu_leaf, sk_leaf)
+
+
+def test_deep_tree_classifier(clf_data):
+    """Deep trees work correctly (depth 10+)."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=12, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    est = rf.estimators_[0]
+    proba = est.predict_proba(X)
+    assert proba.shape == (X.shape[0], 2)
+    np.testing.assert_allclose(proba.sum(axis=1), 1.0, rtol=1e-5)
+    assert est.tree_.max_depth <= 12
+
+
+def test_estimators_cache_invalidated_on_refit(clf_data):
+    """Re-fitting clears the estimators_ cache."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=4, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+    old_est = rf.estimators_
+
+    rf.fit(cp.asarray(X), cp.asarray(y))
+    new_est = rf.estimators_
+    assert old_est is not new_est
+
+
+# ---------------------------------------------------------------------------
+# Ecosystem integration tests
+# ---------------------------------------------------------------------------
+
+
+def test_skforecast_pattern(clf_data):
+    """Mimics skforecast's access: estimators_[i].tree_.predict(X)."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    for i in range(len(rf.estimators_)):
+        tree_pred = rf.estimators_[i].tree_.predict(X[:10])
+        assert tree_pred.shape[0] == 10
+        assert tree_pred.ndim == 2
+        np.testing.assert_allclose(tree_pred.sum(axis=1), 1.0, rtol=1e-5)
+
+
+def test_skforecast_pattern_regressor(reg_data):
+    """Mimics skforecast's access for regression: estimators_[i].tree_.predict(X)."""
+    X, y = reg_data
+    rf = cuRFR(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    for i in range(len(rf.estimators_)):
+        tree_pred = rf.estimators_[i].tree_.predict(X[:10])
+        assert tree_pred.shape[0] == 10
+
+
+def test_decision_path(clf_data):
+    """decision_path returns a sparse matrix with correct shape."""
+    X, y = clf_data
+    rf = cuRFC(n_estimators=3, max_depth=5, random_state=42)
+    rf.fit(cp.asarray(X), cp.asarray(y))
+
+    tree = rf.estimators_[0].tree_
+    path = tree.decision_path(X[:10])
+
+    assert path.shape == (10, tree.node_count)
+    # Each sample must visit the root
+    assert np.all(path[:, 0].toarray() == 1)
+    # Each sample visits at least one leaf
+    is_leaf = tree.children_left == -1
+    for i in range(10):
+        visited = path[i].toarray().ravel()
+        assert np.any(visited[is_leaf])