Как я могу создать образец из логарифмического распределения с хвостом Парето в R?

Question

Я хотел бы создать образец из нормального логарифмического дистрибутива с хвостом Парето в R. Кто-нибудь может мне помочь? Благодарю.

Answer 1

Я не уверен, что это то, что вы ищете, но есть много литературы на тему нормальных распределений с двойным логарифмическим парето или так называемого dPlN. Видишь это, и это, и это. Они используются для моделирования всего: от распределения доходов и доходов до размера нефтяного месторождения и интернет-трафика.

Если это не то, что вы ищете, дайте мне знать, и я удалю сообщение.

Вы спрашиваете, как генерировать случайную выборку из dPlN. Теоретически, можно генерировать случайную выборку из произвольного распределения, генерируя случайную выборку из равномерного распределения, U[0,1] и преобразование этого с использованием обратного CDF целевого распределения.

Итак, сначала нам нужен PDF-файл dPlN, затем мы интегрируем его, чтобы найти CDF, а затем инвертируем его, чтобы найти обратный CDF. PDF dPlN приведен в первой ссылке в уравнениях 8 и 9:

где α и β - параметры местоположения, а ν и τ² - среднее значение и дисперсия логнормального распределения. Φ и Φ^c являются CDF и дополнительным CDF N[0,1]. Грубо говоря, меньшие α и β означают более длинный хвост, ν влияет на положение пика, а τ влияет на ширину пика.

Таким образом, в R мы вычисляем PDF, CDF и обратный CDF dPlN следующим образом:

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, поэтому нет никаких гарантий, что он будет работать при любых обстоятельствах.