在 Python 中模拟相关的对数正态答案

【问题标题】：Simulating correlated lognormals in Python在 Python 中模拟相关的对数正态
【发布时间】：2019-10-30 21:32:25
【问题描述】：

我正在关注这个问题How can I sample a multivariate log-normal distribution in 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，但它有点偏离。由于值已经收敛，因此增加模拟次数没有任何作用。

【问题讨论】：

标签： python numpy random statistics

【解决方案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）。详细讨论了三个变量相关约束can be found on Mathematics 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 几乎可以做它想做的事；幸运的是，它似乎离得很近。

【讨论】：