找到一種方法來模擬這個分佈的隨機數

我正在嘗試在 R 中編寫一個程序，該程序模擬具有累積分佈函數的分佈中的偽隨機數：

$$ F(x)= 1-\exp \left(-ax-\frac{b}{p+1}x^{p+1}\right), \quad x \geq 0 $$

在哪裡 $ a,b>0, p \in (0,1) $

我嘗試了逆變換採樣，但逆變換似乎無法解析解。如果您能提出解決此問題的方法，我會很高興

這個練習有一個簡單的（如果我可以補充的話，優雅的）解決方案：因為 $ 1-F(x) $ 看起來像兩個生存分佈的乘積： $$ (1-F(x))=\exp\left{-ax-\frac{b}{p+1}x^{p+1}\right}=\underbrace{\exp\left{-ax\right}}{1-F_1(x)}\underbrace{\exp\left{-\frac{b}{p+1}x^{p+1}\right}}{1-F_2(x)} $$ 分佈 $ F $ 是分佈 $$ X=\min{X_1,X_2}\qquad X_1\sim F_1,,X_2\sim F_2 $$ 在這種情況下 $ F_1 $ 是指數 $ \mathcal{E}(a) $ 分佈和 $ F_2 $ 是個 $ 1/(p+1) $ - 指數的冪 $ \mathcal{E}(b/(p+1)) $ 分配。

相關的 R 代碼非常簡單

x=pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))) #simulating an n-sample

它絕對比反pdf和accept-reject分辨率快得多：

> n=1e6
> system.time(results <- Vectorize(simulate,"prob")(runif(n)))
utilisateur     système      écoulé 
   89.060       0.072      89.124 
> system.time(x <- simuF(n,1,2,3))
utilisateur     système      écoulé 
    1.080       0.020       1.103 
> system.time(x <- pmin(rexp(n,a),rexp(n,b/(p+1))^(1/(p+1))))
utilisateur     système      écoulé 
    0.160       0.000       0.163 

毫不奇怪的完美契合：

引用自：https://stats.stackexchange.com/questions/391054