  Accurate and precise covariance matrices will be important in enabling
planned cosmological surveys to detect new physics. Standard methods imply
either the need for many N-body simulations in order to obtain an accurate
estimate, or a precise theoretical model. We combine these approaches by
constructing a likelihood function conditioned on simulated and theoretical
covariances, consistently propagating noise from the finite number of
simulations and uncertainty in the theoretical model itself using an
informative Inverse-Wishart prior. Unlike standard methods, our approach allows
the required number of simulations to be less than the number of summary
statistics. We recover the linear 'shrinkage' covariance estimator in the
context of a Bayesian data model, and test our marginal likelihood on simulated
mock power spectrum estimates. We conduct a thorough investigation into the
impact of prior confidence in different choices of covariance models on the
quality of model fits and parameter variances. In a simplified setting we find
that the number of simulations required can be reduced if one is willing to
accept a mild degradation in the quality of model fits, finding that even
weakly informative priors can help to reduce the simulation requirements. We
identify the correlation matrix of the summary statistics as a key quantity
requiring careful modelling. Our approach can be easily generalized to any
covariance model or set of summary statistics, and elucidates the role of
hybrid estimators in cosmological inference.
