6

我是 R 的新手,需要有关如何从具有参数 s 和 beta 的帕累托分布的有限区域中提取随机数的建议。(系统:Windows 7,R 2.15.2。)

(1) 我在向量 data$t 中有数据;每个数据点我都会调用 data&tx

对于这些数据,帕累托分布的参数 s 和 beta 估计如下https://stats.stackexchange.com/questions/27426/how-do-i-fit-a-set-of-data-to-a-pareto -分布在-r

pareto.MLE <- function(X)
{
n <- length(X)
m <- min(X)
a <- n/sum(log(X)-log(m))
return( c(m,a) ) 
}

(2) 现在我需要在这个帕累托分布 (s, beta, 请参见 (1)) 上绘制与观察 (= data points: data$tx) 一样多的随机数 (RndNew) 。对于抽奖,抽取随机数的区域必须限制在 RndNewx >= data$tx; 的区域内。换句话说:RndNewx 绝不能小于相应的 data$tx。

问题:如何告诉 R 限制帕累托分布的区域,从中抽取随机数为 RndNewx >= data$tx?

感谢一百万的任何帮助!

4

1 回答 1

18

从截断分布中抽样的标准方法包括三个步骤。这是正态分布的示例,因此您可以理解。

n <- 1000
lower_bound <- -1
upper_bound <- 1

将 CDF 应用于您的下限和上限,以找到分布边缘的分位数。

(quantiles <- pnorm(c(lower_bound, upper_bound)))
# [1] 0.1586553 0.8413447

从这些分位数之间均匀分布的样本。

uniform_random_numbers <- runif(n, quantiles[1], quantiles[2])

应用逆 CDF。

truncated_normal_random_numbers <- qnorm(uniform_random_numbers)

截断的正态随机数


帕累托分布的 CDF 为

ppareto <- function(x, scale, shape)
{
  ifelse(x > scale, 1 - (scale / x) ^ shape, 0)
}

反之亦然

qpareto <- function(y, scale, shape)
{
  ifelse(
    y >= 0 & y <= 1,
    scale * ((1 - y) ^ (-1 / shape)),
    NaN
  )
}

我们可以修改上面的例子来使用这些帕累托函数。

n <- 1000
scale <- 1
shape <- 1
lower_bound <- 2
upper_bound <- 10

(quantiles <- ppareto(c(lower_bound, upper_bound), scale, shape))
uniform_random_numbers <- runif(n, quantiles[1], quantiles[2])    
truncated_pareto_random_numbers <- qpareto(uniform_random_numbers, scale, shape)

截断帕累托分布


为了更容易重用这段代码,我们可以将它包装到一个函数中。下限和上限具有与分布范围相匹配的默认值,因此如果您不传递值,那么您将获得非截断的帕累托分布。

rpareto <- function(n, scale, shape, lower_bound = scale, upper_bound = Inf)
{
  quantiles <- ppareto(c(lower_bound, upper_bound), scale, shape)
  uniform_random_numbers <- runif(n, quantiles[1], quantiles[2])    
  qpareto(uniform_random_numbers, scale, shape)
}
于 2013-01-24T09:23:56.877 回答