Monte Carlo 使用 Gamma 分布估计 Theta

Question

我想运行在 r 中进行 monte carlo 模拟来估计 theta。有人可以推荐一些资源并建议我如何做到这一点吗？

我已经开始使用 gamma 分布创建样本并使用分布的形状和比率，但我不确定下一步该怎么做。

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)

Answer 1

根据定义估计 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)

不太清楚您希望如何使用 Monte Carlo 模拟来帮助您解决这些问题。

Answer 2

您可以重写 θ 的表达式，分解出 exponential distribution。

θ = ₀∫^∞ (x^4.5/2) (2 e^-2x) dx

这里(2 e^-2x)是rate=2的指数分布，提示如何用Monte Carlo.

积分

来自指数的样本随机值
计算函数 (x^4.5/2) 在采样值
这些计算值的平均值将是 M-C 计算的积分

代码，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) 是我们的幂函数 (x^4.5/2)，f(x) 是我们的指数分布。

更新

只是为了澄清一件事 - 没有单一的规范方法可以将下积分表达式拆分为采样 PDF f(x) 和可计算函数 g(x)，其平均值将是我们的积分。

例如，我可以写

θ = ₀∫^∞ (x^4.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.5^5.5 ₀∫^∞ 1 G(x | shape=5.5, scale=0.5) dx

其中 Γ(x) 是伽马函数，G(x| shape, scale) 是伽马分布。因此，您可以从伽马分布和 g(x)=1 中对任何采样 x 进行采样。因此，这将为您提供准确的答案。代码

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的积分值。

Monte Carlo 使用 Gamma 分布估计 Theta

Monte Carlo to Estimate Theta using Gamma Distribution

simulation

r

gamma

gamma-distribution

montecarlo