我想在 r 中运行蒙特卡罗模拟来估计 theta。有人可以推荐一些资源并建议我如何做到这一点吗？
我已经开始创建具有伽马分布的样本并使用分布的形状和速率，但我不确定下一步该去哪里。
x = seq(0.25, 2.5, by = 0.25)
PHI <- pgamma(x, shape = 5.5, 2)
CDF <- c()
n= 10000
set.seed(12481632)
y = rgamma(n, shape = 5.5, rate = 2)
【问题讨论】：
标签： r simulation montecarlo gamma-distribution gamma
你可以重写你的 θ 表达式，分解出exponential distribution。
θ = 0∫∞ (x4.5/2) (2 e-2x) dx
这里 (2 e-2x) 是 rate=2 的指数分布，这表明如何使用 Monte Carlo 对其进行积分。
代码，R 4.0.3 x64，Win 10
set.seed(312345)
n <- 10000
x <- rexp(n, rate = 2.0)
f <- 0.5*x**4.5
mean(f)
打印
[1] 1.160716
您甚至可以将统计误差估计为
sd(f)/sqrt(n)
打印出来的
[1] 0.1275271
因此，积分 θ 的 M-C 估计为 1.160716∓0.1275271
这里实现的内容如下，例如http://www.math.chalmers.se/Stat/Grundutb/CTH/tms150/1112/MC.pdf，6.1.2，其中 g(x) 是我们的幂函数 (x4.5/2)，f(x) 是我们的指数分布。
更新
只是为了澄清一件事 - 没有单一的规范方法可以将积分表达式拆分为采样 PDF f(x) 和可计算函数 g(x)，它们的平均值将是我们的积分。
例如，我可以写
θ = 0∫∞ (x4.5 e-x) (e-x ) dx
e-x 将是 PDF f(x)。 rate=1 的简单指数，而 g(x) 如何具有指数剩余部分。类似代码
set.seed(312345)
n <- 10000
f <- rexp(n, rate = 1.0)
g <- f**4.5*exp(-f)
print(mean(g))
print(sd(g)/sqrt(n))
产生的整数值为 1.148697∓0.02158325。这是一种更好的方法，因为统计误差更小。
你甚至可以写成
θ = Γ(5.5) 0.55.50∫∞ 1 G(x| shape=5.5, scale=0.5) dx
其中 Γ(x) 是 gamma 函数，G(x| shape, scale) 是 Gamma 分布。因此，您可以从伽马分布中采样，并且对于任何采样的 x，g(x)=1。因此，这将为您提供准确的答案。代码
set.seed(312345)
f <- rgamma(n, 5.5, scale=0.5)
g <- f**0.0 # g is equal to 1 for any value of f
print(mean(g)*gamma(5.5)*0.5**5.5)
print(sd(g)/sqrt(n))
产生的整数值为 1.156623∓0。
【讨论】：
在给定定义的情况下估计 theta 的最佳方法是
theta <- integrate(function(x) x^4.5 * exp(-2*x), from = 0, to = Inf)
给予：
theta
#> [1] 1.156623
另一种处理方法是查看常数 lambda^rate / gamma(rate) 可以在 cdf 积分之外取值，并且由于我们知道无穷远处的 cdf 为 1，那么 theta 必须等于 gamma(rate)/lambda^rate
gamma(5.5)/2^5.5
#> [1] 1.156623
请注意，我们还可以为您的 pdf 和 cdf 编写函数并绘制它们：
pdf <- function(t, rate, lambda) {
(lambda^rate)/gamma(rate) * t^(rate-1) * exp(-2 * t)
}
cdf <- function(x, rate, lambda) {
sapply(x, function(y) {
integrate(pdf, lower = 0, upper = y, lambda = lambda, rate = rate)$value
})
}
curve(pdf(x, 5.5, 2), from = 0, to = 10)
curve(cdf(x, 5.5, 2), from = 0, to = 10)
目前还不清楚您希望蒙特卡洛模拟如何帮助您解决这些问题。
【讨论】：