The effectiveness of variable fiber spacing composite (VFSC) laminates subjected to thermal loading is compared with uniformly distributed fiber composite (UDFC) laminates. The higher-order shear deformation theory (HSDT) is adopted to develop a FE model for estimating critical thermal buckling load. Three types of fiber distribution patterns are considered to study the effect of the non-uniform spacing distribution of fiber on thermal buckling, and the results are compared with UDFC. Further, a surrogate model based on the radial basis function network (RBFN) is developed to examine the impact of composite properties uncertainty on the critical buckling response of UDFC and VFSC laminates under thermal loading. The RBFN model is found to be highly efficient and possesses a close match with Monte Carlo simulation (MCS). Though the VFSC laminates possess a large deviation (standard deviation) in the buckling temperature, the stochastic effect (standard deviation/mean) shows relatively minimal variation in VFSC and UDFC laminates. The effect of each composite property is investigated for different lamination sequences and boundary conditions. From the sensitivity analysis, the thermal expansion coefficient of matrix [Formula: see text] is found to be the most sensitive property for all the investigated cases. Lastly, the RBFN model is utilized to quantify reliability parameters for multiple stochasticity levels.