4

我想从 R 中带有帕累托尾的对数正态分布生成一个样本。有人可以帮我吗？谢谢。

4

1 回答 1

10

我不确定这是否是您要查找的内容，但有大量关于双帕累托对数正态分布或 so-ca.led dPlN 主题的文献。看到这个，这个，这个。这些用于模拟从收入分配、油田规模到互联网流量的一切。

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

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

所以首先，我们需要 dPlN 的 PDF，然后我们整合它以找到 CDF，然后我们反转它以找到逆 CDF。dPlN 的 PDF 在 Eqn 8 和 9的第一个参考文献中给出：

其中 α 和 β 是位置参数，ν 和 τ 2是对数正态分布的均值和方差。Φ 和 Φ 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 的所有可能值进行测试，因此不能保证它在所有情况下都有效。

于 2014-04-22T19:19:25.627 回答