Я хотел бы запустить симуляцию Монте-Карло в r, чтобы оценить тета. Может ли кто-нибудь порекомендовать некоторые ресурсы и предложить, как я могу это сделать?
Я начал с создания выборки с гамма-распределением и с использованием формы и скорости распределения, но я не уверен, что делать дальше.
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)
Вы можете переписать выражение для θ, исключив экспоненциальное распределение.
θ = 0 ∫ ∞ (x 4.5 / 2) (2 e -2x) dx
Здесь (2 e -2x) - экспоненциальное распределение со скоростью = 2, которое подсказывает, как интегрировать его с помощью Монте-Карло.
Код, 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), среднее значение которой было бы нашим интегралом.
Например, я мог бы написать
θ = 0 ∫ ∞ (x 4.5 e -x) (e -x ) dx
e -x будет PDF f (x). Простая экспонента со скоростью = 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 0 ∫ ∞ 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.
Наилучший способ оценить тэту с учетом ее определения -
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, тогда тета должна равно гамма (коэффициент) / лямбда ^ коэффициент
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)
Не совсем понятно, как вы хотите, чтобы моделирование Монте-Карло помогло вам в этом.
