  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 chaotic systems. 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 number
of 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\'{e}nard convection.
