等高线图双变量对数正态密度函数python答案

【问题标题】：contour plot bivariate lognormal density function python等高线图双变量对数正态密度函数python
【发布时间】：2021-02-14 06:18:22
【问题描述】：

我想使用此 python 代码绘制随机变量 R~LN(7, 0.5) 和 S~LN(1, 0.5) 的二元对数正态 PDF 的等高线图：

import numpy as np
import matplotlib.pyplot as plt
from scipy.stats import lognorm

r= np.linspace(1, 10, 500)
s= np.linspace(1, 10, 500)
R, S = np.meshgrid(r, s)

sigma_R = 0.5
mu_R = 7

sigma_S = 0.5
mu_S = 1

#lognormal PDF
pdf_R = lognorm.pdf(R.flatten(), 1, mu_R, sigma_R)
pdf_S = lognorm.pdf(S.flatten(), 1, mu_S, sigma_S)
JointPDF = pdf_R*pdf_S

fig, ax = plt.subplots()
CS2 = ax.contour(R, S, JointPDF.reshape(500,500), 30, cmap="RdBu_r")

结果是：

Q1：这个情节正确吗？我不确定，因为等高线图不应该从 1 和 7 开始，这对应于 RV 的平均值。

Q2：有谁知道万一R和S相关怎么办？

谢谢

【问题讨论】：

标签： python matplotlib

【解决方案1】：

让我们从第一季度开始。不，情节可能不是您想要的，主要是由于 scipy 中 lognorm 的界面有些混乱，并且您的符号 LN(7,0.5) 不清楚。因此，如果您能澄清这将是有帮助的。但为了取得进展，我假设您想要对数正态分布，其均值和 sigma 如 Wikipedia 中所定义，维基百科的参数 mu，sigma 为 mu_R=log(7) 和 sigma_R = 0.5（否则规模的数字没有意义）

所以如果 mu、sigma 是1 中定义的两个参数，那么你要调用lognorm 如下

pdf_R = lognorm.pdf(R.flatten(), sigma_R, 0, mu_R)

（对于 pdf_S 类似）

您还想对您的域小心一点。分布将以 (7,1) 为中心，在 (0,inf) x (0, inf) 上定义，宽度约为 sigma 的几个倍数。如果你想要它是对称的，我会让域 [0,10]x[0,10]:

r= np.linspace(0, 10, 500)
s= np.linspace(0, 10, 500)

Q2 完全是另一回事，所以我将把它留在那个地方

【讨论】：

感谢您的回答。 R 是参数 mu_R=7 和 sigma_R=0.5 的对数正态 RV，所以我应该这样使用 scipy：pdf_R = lognorm.pdf(R.flatten(), sigma_R, 0, exp(mu_R))
这对我来说似乎是对 LN(7,0.5) 的可能解释，唯一的犹豫是 exp(7) 是一个相当大的数字——比 exp(mu_S ) = exp(1) （大约 x400 大）所以你的等高线图看起来真的很奇怪
另一种可能的解释是 LN 变量的预期值为 7，标准差为 0.5（因此方差为 0.5^2=0.25）这将转化为另一组 mu、sigma（使用我链接的维基百科文章中的公式）。
第二种解释的可能性更大。将 mu 和 sigma 的维基百科公式与 scipy 的对数正态一起使用可以得到所需的结果，请注意：pdf_R = lognormal.pdf(R.flatten(), sigma, 0, exp(mu))

【解决方案2】：

作为对我问题 Q2 的回答，我在 https://reference.wolfram.com/language/ref/LogMultinormalDistribution.html 中找到了以下表达式

这是双变量 lognoraml PDF 的显式公式。这是一个初学者代码，但它可能会有所帮助

import numpy as np
import matplotlib.pyplot as plt

""" lognormal parameters """
def mu_z(mu, sigma):
    mu_Z = np.log(mu/(np.sqrt(1+(sigma/mu)**2)))
    return mu_Z
def sigma_z(mu, sigma):
    sigma_Z = np.sqrt(np.log(1+(sigma/mu)**2))
    return sigma_Z

r= np.linspace(1, 10, 500)
s= np.linspace(1, 10, 500)
R, S = np.meshgrid(r, s)

sigma_R = 0.5
mu_R = 7
Z_R = np.log(R)
sigma_ZR = sigma_z(mu_R, sigma_R)
mu_ZR = mu_z(mu_R, sigma_R)

sigma_S = 0.5
mu_S = 1
Z_S = np.log(S)
sigma_ZS = sigma_z(mu_S, sigma_S)
mu_ZS = mu_z(mu_S, sigma_S)

rho = 0.5
A = 1/(2 * R.flatten()* S.flatten()*np.pi*np.sqrt(sigma_ZR**2*sigma_ZS**2-rho**2*sigma_ZR**2*sigma_ZS**2))
C = (np.log(R.flatten())-mu_ZR)*(-((rho*(np.log(S.flatten())-mu_ZS)*sigma_ZR*sigma_ZS)/
    (sigma_ZR**2*sigma_ZS**2-rho**2*sigma_ZR**2*sigma_ZS**2))+
    (((np.log(R.flatten())-mu_ZR)*sigma_ZS**2)/
    (sigma_ZR**2*sigma_ZS**2-rho**2*sigma_ZR**2*sigma_ZS**2)))
B = (np.log(S.flatten())-mu_ZS)*((((np.log(S.flatten())-mu_ZS)*sigma_ZR**2)/
    (sigma_ZR**2*sigma_ZS**2-rho**2*sigma_ZR**2*sigma_ZS**2))-
    ((rho*(np.log(R.flatten())-mu_ZR)*sigma_ZS*sigma_ZR)/
    (sigma_ZR**2*sigma_ZS**2-rho**2*sigma_ZR**2*sigma_ZS**2)))
    
pdflnmv = A*np.exp(-(C+B)/2)
fig, ax = plt.subplots()
CS4 = ax.contour(R, S, pdflnmv.reshape(500,500))

当rho = 0.5（rho为相关系数）时，结果为

当 rho = 0.0（不相关的 RV）时，我们就有了 Q1 的答案

【讨论】：