We present a flexible data-driven framework for estimating the response of higher-order moments of nonlinear stochastic systems to small external perturbations. The classical generalized fluctuation–dissipation theorem (GFDT) links the unperturbed steady-state distribution to the system’s linear response. While standard implementations relying on Gaussian approximations can predict the mean response, they often fail to capture changes in higher-order moments. To overcome this, we combine GFDT with score-based generative modeling to estimate the system’s score function directly from data. We demonstrate the framework’s versatility by employing two complementary score estimation techniques tailored to the system’s characteristics: i) a clustering-based algorithm (K-means Gaussian Mixture Modeling) for systems with low-dimensional effective dynamics, and ii) a denoising score matching method implemented with a U-Net architecture for high-dimensional, spatially extended systems where reduced-order modeling is not feasible. Our method is validated on several stochastic models relevant to climate dynamics: three reduced-order models of increasing complexity and a 2D Navier–Stokes model representing a turbulent flow with a localized perturbation. In all cases, the approach accurately captures strongly nonlinear and non-Gaussian features of the system’s response, significantly outperforming traditional Gaussian approximations.