Data analysis in cosmology requires reliable covariance matrices. Covariance matrices derived from numerical simulations often require a very large number of realizations to be accurate. When a theoretical model for the covariance matrix exists, the parameters of the model can often be fit with many fewer simulations. We establish a rigorous Bayesian method for performing such a fit, but show that using the maximum posterior point is often sufficient. We show how a model covariance matrix can be tested by examining the appropriate χ2 distributions from simulations. We demonstrate our method on two examples. First, we measure the two-point correlation function of halos from a large set of 10000 mock halo catalogs. We build a model covariance with 2 free parameters, which we fit using our procedure. The resulting best-fit model covariance obtained from just 100 simulation realizations proves to be as reliable as the numerical covariance matrix built from the full 10000 set. We also test our method on a setup where the covariance matrix is large by measuring the halo bispectrum for thousands of triangles for the same set of mocks. We build a block diagonal model covariance with 2 free parameters as an improvement over the diagonal Gaussian covariance. Our model covariance passes the χ2 test only partially in this case, signaling that the model is insufficient even using free parameters, but significantly improves over the Gaussian one.