This paper presents a new probabilistic method based on Hybrid High Dimensional Model Representation (HHDMR) for predicting the failure probability of structural and mechanical systems subject to random loads, material properties, and geometry. High Dimensional Model Representation (HDMR) is a general set of quantitative model assessment and analysis tools for capturing the high dimensional relationships between sets of input and output model variables. It is a very efficient formulation of the system response, if higher order variable cooperative effects are weak and if the response function is dominantly of additive nature, allowing the physical model to be captured by the first few lower order terms. But, if multiplicative nature of the response function is dominant then all right hand side components of HDMR must be used to be able to obtain the best result. However, if HDMR requires all components, which means 2 N number of components, to get a desired accuracy, making the method very expensive in practice, then Factorized HDMR (FHDMR) can be used. But in most cases the limit state/performance function has neither additive nor multiplicative nature, rather it has an intermediate nature. This paper presents a new HHDMR-based approximation for the limit state/performance functions having an intermediate nature. The proposed approximation of an implicit limit state/performance function includes both HDMR and FHDMR expansions through a hybridity parameter. As an alternative to the conventional methods for reliability analysis which are very computationally demanding, when applied in conjunction with complex finite element models, this study aims to assess how accurately and efficiently HHDMR technique can capture complex model output uncertainty. Once the approximate form of the original implicit limit state/performance function is defined, the failure probability can be obtained by statistical simulation. Results of six numerical examples involving mathematical functions and structural/solid-mechanics problems indicate that the failure probability obtained using HHDMR-based approximation of an implicit limit state/performance function provides significant accuracy when compared with the conventional Monte Carlo method, while requiring fewer original model simulations.