如何找到一组数据的分布，然后进一步传播此分布

我的问题是，我有一组数据，我想将分布拟合到其中，然后在找到分布后，对其运行蒙特卡罗模拟，以传播找到的分布。

我的第一个代码是：

require(fitdistrplus)
example1<-c(29,23,29,25,26,29,29,27,25,25,25,26,28,25,29,28,28,26,28,25,29,26,30) 
f1<-fitdist(example1,rgamma,method="mle")

如果我使用命令

print(f1)

它告诉我形状是204.00伽马分布的速率为7.568

（请注意，目前我拟合分布的数字是任意的，我通常有数百个观测值来拟合分布）。

我现在需要帮助的地方是，当我使用包mc2d中的代码来传播此分发时，如下所示：

require(mc2d)
ndunc(1000)
fitted<-mcstoc(rgamma, type="U", shape=204.00, rate=7.569)

目前，我不得不从"fitdist"命令的前一个"打印"中手动输入形状和速率到上面的函数中。

我的问题是，有没有一种方法可以让mcstoc命令自动从fitdist命令中获取形状和速率，这样我就不必手动中断代码了？或者，如果fitdistrplus软件包和mc2d软件包不可能，那么还有其他软件包可以为我做到这一点吗？

非常感谢！

f1$estimate[1]
#   shape 
#204.0008 
f1$estimate[2]
#    rate 
#7.567762
fitted<-mcstoc(rgamma, type="U", shape=f1$estimate[1], rate=f1$estimate[2])

myFunction <- function (data){
    f1<-fitdist(data,rgamma,method="mle")
    fitted<-mcstoc(rgamma, type="U", shape=f1$estimate[1], rate=f1$estimate[2])
    return(fitted)
}
example1<-c(29,23,29,25,26,29,29,27,25,25,25,26,28,25,29,28,28,26,28,25,29,26,30)
fitted.example1 <- myFunction(exemple1)

fitted <- do.call( 
  function(...) mcstoc(rgamma, type="U", ...), 
  as.list(f1$estimate) 
)