Имитация коррелированных логнормалей в Python

Question

Я следую за ответом в этом вопросе Как я могу сэмплировать многомерное логарифмическое нормальное распределение в Python? , но я получаю, что предельные распределения выборочных данных не могут иметь одинаковыесреднее и стандартное отклонение введенных маргиналов. Например, рассмотрим многомерное распределение ниже в примере кода. Если мы помечаем маргиналы как X, Y и Z, то я ожидаю, что параметры масштаба и местоположения (подразумеваемые из данных выборки) будут соответствовать введенным данным. Однако для X вы можете видеть ниже, что параметры масштаба и местоположения равны 0,1000 и 0,5219. Таким образом, масштаб - это то, что мы ожидаем, но местоположение отключено на 4%. Я думаю, что я делаю что-то не так с ковариационной матрицей, но я не могу понять, что не так. Я попытался установить матрицу корреляции в единичную матрицу, а затем расположение и масштаб данных образца совпадают с введенными данными. Что-то должно быть не так с моей ковариационной матрицей, или я делаю еще одну фундаментальную ошибку. Любая помощь будет оценена. Пожалуйста, сообщите, если вопрос неясен.

import pandas as pd
import numpy as np
from copy import deepcopy

mu  = [0.1, 0.2, 0.3]
sigma = [0.5, 0.8, 0.6]
sims = 3000000
rho = [[1, 0.9, 0.3], [0.9, 1, 0.8], [0.3, 0.8 ,1]]

cov = deepcopy(rho)
for row in range(len(rho)):
    for col in range(len(rho)):
        cov[row][col] = rho[row][col] * sigma[row] * sigma[col]

mvn = np.random.multivariate_normal(mu, cov, size=sims) 

sim = pd.DataFrame(np.exp(mvn), columns=['X', 'Y', 'Z'])

def computeImpliedLogNormalsParams(mean, std):
    # This method implies lognormal params which match the moments inputed 
    secondMoment = std * std + mean *mean
    location = np.log(mean*mean / np.sqrt(secondMoment))
    scale = np.sqrt(np.log(secondMoment / (mean * mean)))
    return (location, scale)

def printDistributionProp(col, sim):
    print(f"Mean = {sim[col].mean()}, std = {sim[col].std()}")
    location, scale = computeImpliedLogNormalsParams(sim[col].mean(), sim[col].std())
    print(f"Matching moments gives a lognormal with location {location} and scale {scale}")


printDistributionProp('X', sim)

Вывод:

Mean = 1.2665338803521895, std = 0.708713940557892
Matching moments gives a lognormal with location 0.10008162992913544 and scale 0.5219239625443672   

Наблюдая за выводом, мы ожидаем, что параметр масштаба будет очень близок к 0,5, но он немного отключен. Увеличение количества симуляций ничего не дает, так как значение сблизилось.

Answer 1

Ковариационная матрица не является положительной полуопределенной:

>>> mvn = np.random.multivariate_normal(mu, cov, size=sims, check='raise')
Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "mtrand.pyx", line 4542, in mtrand.RandomState.multivariate_normal
ValueError: covariance is not symmetric positive-semidefinite.

и, следовательно, не существует распределения данных, которое фактически имеет запрошенную ковариационную структуру. На высоком уровне учтите, что вы указываете X и Z , чтобы оба были сильно коррелированы с Y (0,8 и 0,9), но в то же времябыть довольно слабо коррелированными друг с другом (0,3). Подробное обсуждение, в частности, трех переменных корреляционных ограничений , можно найти в Математике SE .

Я не знаю, как NumPy справляется с этим (вы должны были увидеть предупреждение),но если вы проверите окончательную структуру корреляции:

>>> np.corrcoef(mvn.T)
array([[1.        , 0.79817321, 0.33343102],
       [0.79817321, 1.        , 0.74525583],
       [0.33343102, 0.74525583, 1.        ]])

, можно увидеть, что X и Z имеют более низкие корреляции с Y иболее высокая корреляция друг с другом, чем изначально указано rho. Опять же, не уверен, как именно корректируются отклонения, но поскольку ковариация невозможна, NumPy может в значительной степени делать то, что хочет;к счастью, кажется, что он довольно близко.