Machine learning (ML) has emerged as a powerful approach for accelerating the discovery of thermoelectric (TE) materials, yet many models suffer from limited generalizability. Here we present a robust ML framework for predicting the thermoelectric figure of merit (zT) in half-Heusler (hH) compounds, emphasizing prediction of previously unseen compositions. A principal component analysis (PCA)-based train/test splitting strategy is introduced to generate representative partitions of chemical space. Bayesian hyperparameter optimization is combined with correlation-based k-best feature selection for Random Forest, XGBoost, and Neural Network models, together with the interpretable SISSO symbolic-regression method. Relative to the full set of 118 input features, the optimized feature subsets reduce the test RMSE by 16-27%. SHAP and SISSO analyses identify temperature, A-site dopant concentration, and A-site heat of vaporization as the most important descriptors of zT. The framework is validated through prediction of previously unexplored compositions, recovery of literature-reported optimal zT values for established hH systems, and first-principles calculations for the new compounds CeCoBi and CeIrBi, which exhibit favorable electronic transport properties. The framework is further applied to high-throughput screening of the doped hH chemical space, covering about 2.46 million elemental substitution combinations across 631 thermodynamically screened parent matrices under stability constraints. This screening yields a concise set of promising candidates for future study. Overall, the work establishes a generalizable workflow for reliable data-driven thermoelectric materials discovery.