Fatigue life estimation, reliability and durability are of paramount importance in acquisition, maintenance and operation of vehicle systems. Fatigue life is random because of the stochastic load, the inherent variability of material properties, and the uncertainty in the definition of the S-N curve. The commonly used fatigue life estimation methods calculate the mean (not the distribution of) fatigue life under Gaussian loads using the potentially restrictive narrow-band assumption. In this paper, we present a general methodology to calculate the distribution of fatigue life of a linear vibratory system under stationary, non-Gaussian loads considering the effects of skewness and kurtosis. The input loads are first characterized using their first four moments (mean, standard deviation, skewness and kurtosis) and a correlation structure equivalent to a given Power Spectral Density (PSD). Then, the first four moments and the correlation structure of the response process are calculated using the impulse response function of the linear vibratory system, and a combination of Polynomial Chaos Expansion (PCE) and Karhunen Loeve (KL) expansion. Finally, the characterized response process is used to estimate the fatigue life using the Miner’s linear damage model and rainflow counting. The approach can be easily extended to non-stationary loads, and with some modifications, to even non-linear vibratory systems. It can also be used to develop an accelerated testing procedure without increasing the applied load artificially. The proposed methodology is demonstrated with an example of a quarter car under non-Gaussian load.