We describe a simple and systematic method for obtaining approximate sensitivity information from a chaotic dynamical system using a hierarchy of cumulant equations. The resulting forward and adjoint systems yield information about gradients of functionals of the system and do not suffer from the convergence issues that are associated with the tangent linear representation of the original chaotic system. The functionals on which we focus are ensemble-averaged quantities, whose dynamics are not necessarily chaotic; hence we analyse the system’s statistical state dynamics, rather than individual trajectories. The approach is designed for extracting parameter sensitivity information from the detailed statistics that can be obtained from direct numerical simulation or experiments. We advocate a data-driven approach that incorporates observations of a system’s cumulants to determine an optimal closure for a hierarchy of cumulants that does not require the specification of model parameters. Whilst the sensitivity information from the resulting surrogate model is approximate, the approach is designed to be used in the analysis of turbulence, whose degrees of freedom and complexity currently prohibits the use of more accurate techniques. Here we apply the method to obtain functional gradients from low-dimensional representations of Rayleigh-Bénard convection.