1

R での再帰は初めてです。R で特定の条件を満たす指数確率変数を生成しようとしています。2 つの独立した指数確率変数を生成する試みを示す簡単な例を次に示します。

コード

#### Function to generate two independent exponential random variable that meet two criteria
gen.exp<-function(q1,q2){
  a=rexp(1,q1)  # generate exponential random variable
  b=rexp(1,q2)  # generate exponential random variable
  if((a>10 & a<10.2) & (a+b>15)){ # criteria the random variables must meet
   return(c(a,b))
  }else{
   return(gen.exp(q1,q2)) #if above criteria is not met, repeat the process again
  } 

}

例: q1=.25、q2=2

     gen.exp(.25,2)

上記のコードを実行すると、次のエラーが発生します。

  1. エラー: 評価の入れ子が深すぎます: 無限再帰 / options(expressions=)?
  2. ラップアップ中のエラー: 評価の入れ子が深すぎます: 無限再帰 / options(expressions=)?

options(expressions=10000)Rが反復回数を増やすことができるようにするために、私はすでに修正を試み ました。それは私の場合には役に立たないようです（オプションを正しく使用していない可能性があります）。上記のような厳しい基準で連続分布を生成することが問題である可能性があることを理解しています。そうは言っても、エラーを回避する方法はありますか? または、少なくともエラーが発生するたびに再帰を繰り返しますか? ここで再帰はオーバーキルですか？目的の確率変数を生成するためのより簡単で優れた方法はありますか? 洞察に感謝します。

4

2 に答える 2

4

2 つの条件で単純な拒否サンプラーを使用している

  • a > 10 & a < 10.2
  • a+ b > 15

ただし、両方が一致する可能性は低く、非常に遅いです。ただし、指数乱数に関心があるため、拒否する数値のシミュレーションを避けることができます。

指数乱数を生成するには、次の式を使用します

-rate * log(U)

は乱数ですU。U(0,1)したがって、指数分布から 10 より大きい (たとえば) 値を生成するには、次のようにします。

-log(U(0, exp(-10*rate))/rate

またはRコードで

-log(runif(1, 0, exp(-10*rate)))/rate

上限にも同様のトリックを使用できます。

上記の @Roland の関数を使用すると、次のようになります。

gen.exp = function(q1, q2, maxiter = 1e3){
  i = 0
  repeat {
    i = i + 1
    upper = exp(-10*q1)
    lower = exp(-10.2*q1)
    a = -log(runif(1, lower, upper))/q1
    b = -log(runif(1, 0, exp(-4.8*q2)))/q2

    if((a>10 & a<10.2) & (a+b>15)) {message(i); return(c(a,b)) 
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries."))
  }
}

繰り返しにかかった回数も出力したことに注意してください。パラメータについては、約2回の反復が必要です。

于 2014-11-27T14:52:07.817 に答える
2

これには再帰を使用しないでください。

gen.exp<-function(q1, q2, maxiter = 1e3){
  i <- 0
  repeat {
    i <- i + 1
    a=rexp(1,q1)  # generate exponential random variable
    b=rexp(1,q2)  # generate exponential random variable
    if((a>10 & a<10.2) & (a+b>15)) return(c(a,b)) # criteria the random variables must meet
    if (i > maxiter) stop(paste("Conditions not fulfilled after", maxiter, "tries."))
  }
}

set.seed(42)
gen.exp(.25, 2)
#Error in gen.exp(0.25, 2) : Conditions not fulfilled after 1000 tries.

bが 4.8 より大きい確率は次のとおりです。

pexp(4.8, 2, lower.tail = FALSE)
#[1] 6.772874e-05

さらに反復してみましょう。

gen.exp(.25, 2, maxiter = 1e7)
#[1] 10.08664  5.55414

もちろん、この RNG は非常に遅いため、ほとんど役に立ちません。一度に大量のaとを生成する方がよいでしょう。b

于 2014-11-27T14:25:39.867 に答える