Неточные результаты при восстановлении ковариационной матрицы из np.random.multivariate_normal - PullRequest
Новая
       43

Неточные результаты при восстановлении ковариационной матрицы из np.random.multivariate_normal

0 голосов
/ 09 января 2019

Мне нужно смоделировать данные из 2-мерного нормального распределения вместе с параметром корреляции. Для этого я использовал np.random.multivariate_normal с ковариационной матрицей, в которой мои квадратные сигмы представлены в виде диагональных элементов и произведений сигм и коэффициента корреляции в других местах (я надеюсь, что это правильный способ получения данных с корреляцией) .

Но, боюсь, я не понимаю, как правильно восстановить ковариационную матрицу из сгенерированных данных. Я пытался получить ковариационную матрицу с np.cov и пытался привести сгенерированные данные к форме с нулевым средним, а затем создать ковариационную матрицу с помощью точечного произведения этих данных.

Вот мой код: 

import numpy as np
from matplotlib import pyplot as plt


class NormalDist:
    def __init__(self, *args):
        self.mu = args[:2]
        self.sigma = args[2:4]
        self.dist, self.cov = None, None

    def generate(self, rho=0., n=100):
        """ generate distributed data """
        self.cov = np.diag(np.array(self.sigma, np.float))
        self.cov = np.power(self.cov, 2)
        corr = rho * self.sigma[0] * self.sigma[1]
        self.cov[0, 1], self.cov[1, 0] = corr, corr
        self.dist = np.random.multivariate_normal(self.mu, self.cov, n)


if __name__ == '__main__':
    gauss = NormalDist(1, 2, 4, 9)
    gauss.generate(1/3)

    # covariance matrix from np.cov
    print(np.cov(gauss.dist.T), '\n')

    # covariance matrix from reducing data to zero-mean form
    zero_mean = gauss.dist - gauss.dist.mean(axis=0, keepdims=True)
    print(zero_mean.T @ zero_mean)

Выход: 

[[13.84078951  9.60607718]
 [ 9.60607718 79.33658308]] 

[[1370.23816181  951.00164066]
 [ 951.00164066 7854.32172506]]

1 Ответ

0 голосов
/ 09 января 2019

вам просто нужно разделить на размер выборки, то есть: 

def np_mv_cov(X):
    X = X - X.mean(axis=0, keepdims=True)
    return (X.T @ X) / (X.shape[0] - 1)

можно протестировать с помощью упрощенной версии кода, указанного выше: 

import numpy as np

dist = np.random.multivariate_normal([1, 2], [[16, 12], [12, 81]], 100)

d = np.cov(dist.T) - np_mv_cov(dist)

print(np.max(np.abs(d)))

дает мне ~ 1.42e-14.

...