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 write a likelihood-based method for performing such a fit. We demonstrate how a model covariance matrix can be tested by examining the appropriate χ 2 distributions from simulations. We show that if model covariance has amplitude freedom, the expectation value of second moment of χ 2 distribution with a wrong covariance matrix will always be larger than one using the true covariance matrix. By combining these steps together, we provide a way of producing reliable covariances without ever requiring running a large number of 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.