r - 如何在 R 中从具有 Pareto 尾部的对数正态分布生成样本？

Question

我想从 R 中具有 Pareto 尾部的对数正态分布生成样本。有人可以帮助我吗？谢谢。

Answer 1

我不确定这就是您要找的内容，但是有大量关于双帕累托对数正态分布或 so-ca.led dPlN 主题的文献。请参阅this ，和this ，和this 。这些用于模拟从 yield 和收入分配到油田规模，再到internet traffic的一切。 .

如果这不是您正在寻找的内容，请告诉我，我将删除该帖子。

您询问如何从 dPlN 生成随机样本。理论上，可以通过从均匀分布 U[0,1] 生成随机样本并使用目标分布的逆 CDF 对其进行转换，从任意分布生成随机样本。

因此，首先，我们需要 dPlN 的 PDF，然后对其进行积分以找到 CDF，然后将其求逆以找到逆 CDF。 dPlN 的 PDF 在 first reference 中给出在等式8和9中:

其中 α 和 β 是位置参数，ν 和 τ² 是对数正态分布的均值和方差。 Φ 和 Φ^c 是 N[0,1] 的 CDF 和互补 CDF。粗略地说，较小的 α 和 β 意味着尾部较长，ν 影响峰的位置，τ 影响峰的宽度。

因此，在 R 中，我们计算 dPlN 的 PDF、CDF 和逆 CDF，如下所示:

f = function(x,alpha, beta, nu, tau) {   # probability density of dPlN
  A = function(theta, nu, tau) exp(theta*nu +(alpha*tau)^2/2)
  c = alpha*beta/(alpha+beta)
  z.alpha = (log(x) - nu - alpha*tau^2)/tau
  z.beta  = (log(x) - nu + beta*tau^2)/tau
  t.alpha = x^-(alpha+1)*A(alpha,nu,tau)*pnorm(z.alpha)
  t.beta  = x^(beta-1)*A(-beta,nu,tau)*(1-pnorm(z.beta))
  return(c*(t.alpha + t.beta))
}
F = function(x,alpha,beta,nu,tau) {      # cumulative density function of dPlN
  ifelse(x > 1e4, 1, integrate(f,0.001,x,alpha,beta,nu,tau)$value)}
F = Vectorize(F, vectorize.args="x")

F.inv = function(y, alpha,beta,nu,tau){  # inverse CDF of dPlN
  uniroot(function(x, alpha,beta,nu,tau){F(x, alpha,beta,nu,tau)-y},
          interval=c(0,1e6),alpha,beta,nu,tau)$root
}
F.inv = Vectorize(F.inv, vectorize.args="y")

x=seq(0,50,length.out=1000)
y=seq(0,.995,length.out=1000)

par(mfrow=c(1,3))
plot(x,f(x,2,2,2,1),type="l",main="f(x)")
plot(x,F(x,2,2,2,1),type="l",main="CDF of f(x)")
plot(y,F.inv(y,2,2,2,1),type="l",main="Inverse CDF of f(x)")

最后，我们使用 F.inv(...) 生成 dPlN 的随机变量，并绘制结果以证明随机样本确实遵循预期的概率分布。

# random sample from dPlN (double Pareto Lognormal distribution)
X = runif(1000,0,1)   # random sample from U[0,1]
Z = F.inv(X,2,2,2,1)

par(mfrow=c(1,1))
hist(Z, breaks=c(seq(min(x),max(x),length=50),Inf), 
     xlim=range(x), freq=FALSE)
lines(x,f(x,2,2,2,1),main="Density function",
      xlim=range(x), col="red", lty=2)

免责声明此代码尚未使用所有可能的 alpha、beta、nu 和 tau 值进行测试，因此不能保证它在所有情况下都能工作。