2013-06-13 4 views
9

私は、同じ分布から来たと思われる一連の数値を持っているとしましょう。密度オブジェクトから乱数を生成する(またはより広範に番号の集合から)

set.seed(20130613) 
x <- rcauchy(10) 

同じ不明な分布からランダムに数値を生成する関数が必要です。私が考えている1つのアプローチは、densityオブジェクトを作成し、そこからCDFを取得し、ランダム一様変数(see Wikipedia)の逆CDFを取ることです。

den <- density(x) 

#' Generate n random numbers from density() object 
#' 
#' @param n The total random numbers to generate 
#' @param den The density object from which to generate random numbers 
rden <- function(n, den) 
{ 
     diffs <- diff(den$x) 
     # Making sure we have equal increments 
     stopifnot(all(abs(diff(den$x) - mean(diff(den$x))) < 1e-9)) 
     total <- sum(den$y) 
     den$y <- den$y/total 
     ydistr <- cumsum(den$y) 
     yunif <- runif(n) 
     indices <- sapply(yunif, function(y) min(which(ydistr > y))) 
     x <- den$x[indices] 

     return(x) 
} 

rden(1, den) 
## [1] -0.1854121 

私の質問は次のとおりです。

  1. は密度オブジェクトから乱数を生成するためのより良い(またはRに組み込まれた)方法はありますか?
  2. 数字のセットから乱数を生成する方法については他にもありますか(sample以外)?
+0

この背後にある理論ははるかに微妙です。密度はどのように推定されますか?どのカーネルが使用されていますか?この見積もりの​​周りに信頼帯がありますか?それは混合モデルですか?等 –

答えて

9

密度推定値からデータを生成するには、元のデータポイントの1つをランダムに選択し、密度推定値からカーネルに基づいてランダムな「エラー」ピースを追加します。デフォルトの「ガウス」では、ランダム元のベクトルからの要素と使用帯域幅に等しい平均0及びSDを有するランダム正常を追加:

den <- density(x) 

N <- 1000 
newx <- sample(x, N, replace=TRUE) + rnorm(N, 0, den$bw) 

別のオプションは、logsplineパッケージからlogspline関数を使用して密度を適合することである(別の方法を使用し密度を推定する)、そのパッケージ内のrlogspline関数を使用して、推定密度から新しいデータを生成します。

2

既存の数値プールから値を引き出すことが必要な場合は、sampleが最適です。
仮定している分布から描画したい場合は、densityを使用し、それを必要な係数(平均、sdなど)を得るために推定分布に適合させ、適切なR分布関数を使用します。

これ以外にも、配布物に応じて「選択的に」サンプリングする方法については、Cの数値レシピの第7.3章(「拒否方法」)を参照してください。コードは簡単にRに簡単に翻訳できるほど簡単です。 私の賭けはすでに誰かがこれを行っているので、これよりも良い答えを投稿します。

関連する問題