We present a scheme for investigating arbitrary thermal observables in spatially inhomogeneous equilibrium many-body systems. Extending the grand canonical ensemble yields any given observable as an explicit hyperdensity functional. Associated local fluctuation profiles follow from an exact hyper-Ornstein-Zernike equation. While the local compressibility and simple observables permit analytic treatment, complex order parameters are accessible via simulation-based supervised machine learning of neural hyperdirect correlation functionals. We exemplify efficient and accurate neural predictions for the cluster statistics of hard rods, square-well rods, and hard spheres. The theory allows one to treat complex observables, as is impossible in standard density functional theory.