时间戳

首页 > 解决方案 > 在 Python 中模拟相关的对数正态

python - 在 Python 中模拟相关的对数正态

问题描述

我正在关注这个问题的答案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，但它有点偏离。由于值已经收敛，因此增加模拟次数没有任何作用。

标签： pythonnumpyrandomstatistics

解决方案

协方差矩阵不是半正定的：

>>> 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 几乎可以做它想做的事；幸运的是，它似乎离得很近。

推荐阅读