我假设您正在参考逆变换采样方法。它非常直截了当。请参阅Wiki文章和此站点。

帕累托 CDF 由下式给出：

$$F(X) = 1 -(\frac{k}{x})^\gamma; x\ge k>0 \ and \ \gamma>0$$

您所做的就是等同于均匀分布并求解$x$

$$F(X) = U \\ U \sim Uniform(0,1) \\ 1 -(\frac{k}{x})^\gamma = u \\ x = k(1-u)^{-1/\gamma}$$

现在在 R 中：

#N = number of samples
#N = number of sample
rpar <- function(N,g,k){
  
  if (k < 0 | g <0){
    stop("both k and g >0")
  }
  
 k*(1-runif(N))^(-1/g)
}


rand_pareto <- rpar(1e5,5, 16)
hist(rand_pareto, 100, freq = FALSE)

#verify using built in random variate rpareto in package extrDistr
x <- (extraDistr::rpareto(1e5,5, 16))
hist(x, 100, freq = FALSE)

这将为您提供帕累托分布的随机变量。我不确定你在哪里得到gammaroot？

使用分位数并反转 CDF 方程给出分位数函数：$p=F(x)$

$$Q(p) = \frac{k}{(1-p)^{1/\gamma}} \quad \quad \quad \text{for all } 0 \leqslant p \leqslant 1.$$

对应的对数分位数函数为：

$$\log Q(p) = \log k - \frac{1}{\gamma} \log (1-p) \quad \quad \quad \text{for all } 0 \leqslant p \leqslant 1.$$

帕累托分布的概率函数已经在R（参见例如EnvStats包）中可用。然而，如果愿意的话，从头开始编程这个函数是相当简单的。这是该函数的矢量化版本。

qpareto <- function(p, scale, shape = 1, lower.tail = TRUE, log.p = FALSE) {
  
  #Check input p
  if (!is.vector(p))                { stop('Error: p must be a numeric vector') }
  if (!is.numeric(p))               { stop('Error: p must be a numeric vector') }
  if (min(p) < 0)                   { stop('Error: p must be between zero and one') }
  if (max(p) > 1)                   { stop('Error: p must be between zero and one') }
  
  n   <- length(p)
  OUT <- numeric(n)
  for (i in 1:n) { 
    P      <- ifelse(lower.tail, 1-p[i], p[i])
    OUT[i] <- log(scale) - log(P)/shape) }
  
   if (log.p) OUT else exp(OUT) }

一旦有了 quantile 函数，就很容易使用逆变换采样生成随机变量。同样，这已经在 中的现有 Pareto 函数中完成R，但如果您想从头开始编程，这很简单。

rpareto <- function(n, scale, shape = 1) {
  
  #Check input n
  if (!is.vector(n))                { stop('Error: n must be a single positive integer') }
  if (!is.numeric(n))               { stop('Error: n must be a single positive integer') }
  if (length(n) != 1)               { stop('Error: n must be a single positive integer') }
  if (as.integer(n) != n)           { stop('Error: n must be a single positive integer') }
  if (n <= 0)                       { stop('Error: n must be a single positive integer') }
  
  qpareto(runif(n), scale = scale, shape = shape) }