The cross-wind response generally follows an un-skewed hardening non-Gaussian distribution around vortex-resonance wind speed, a great amount of samples are required to ensure the accuracy of the peak factor. This study aims to propose a robust methodology for estimating the peak factor, even when dealing with limited samples. The accurate description of the probability distribution for cross-wind response is essential to achieve a reliable peak factor. Since both the harmonic self-excited and Gaussian buffeting components contribute to the non-Gaussianity of the cross-wind response, this study explicates the probability density function (PDF) of the composite process combining the harmonic and Gaussian elements. Subsequently, the PDF of displacement is ascertained by the energy ratio between self-excited vibration and random buffeting, also a kurtosis-based function that can be derived from limited samples. Then the cumulative distribution function (CDF) of the extreme is derived based on the displacement PDF, and a closed-form solution determined by kurtosis for peak factor is derived. Finally, the validity of the proposed method is verified through both Monte Carlo simulations and field measurements. This method offers a practical mean to estimate the extreme values of hardening non-Gaussian responses in highly flexible structures using a limited number of samples.