Ну оценки для элементов ковариационной матрицы -- это просто соответствующие выборочные моменты второго порядка, т.е. просто средние значения соответствующих произведений по всем элементам выборки. Выборочные средние -- тем более. Вот это и надо подставлять в стандартную формулу для плотности многомерного нормального распределения.
Только для несмещённости выборочной ковариационной матрицы делить при усреднении надо не на

, а на

(для матожиданий -- на

).