6

R の初心者で、パラメーター s と beta を使用して、パレート分布の限られた領域から乱数を引き出す方法についてアドバイスが必要です。(システム: Windows 7、R 2.15.2。)

(1) ベクトル data$t にデータがあります。個々のデータ ポイントを data&tx と呼びます

これらのデータについて、パレート分布のパラメーター s とベータは、 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) ここで、観測 (= データ ポイント: data$tx) と同じ数の乱数 (RndNew) をこのパレート分布 (s、ベータ、(1) を参照) から引き出す必要があります。描画の場合、乱数が描画される領域は、RndNewx >= data$tx; の領域に制限する必要があります。つまり、RndNewx は、対応する data$tx より小さくなってはなりません。

質問: RndNewx >= data$tx になるように乱数を引き出すパレート分布の領域を制限するように R に指示するにはどうすればよいですか?

助けてくれてありがとう！

4

1 に答える 1

18

切り捨てられた分布からサンプリングする標準的な方法には、3 つのステップがあります。アイデアを得ることができるように、正規分布の例を次に示します。

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 に答える