Проблемы численного интегрирования через R

Question

У меня есть следующая функция

f(x)∝|x| exp(-1/2 |x| )+1/(1+(x-40)^4 ),xϵR

Я хочу выяснить E(X) и E(X^3) с помощью метода Симпсона (численное интегрирование), стандартного подхода Монте-Карло, выборки с акцептованным отклонением, выборки по важности, алгоритма Метрополиса-Хастинга, выборки Гиббса, а затем байесовской модели с использованием MCMC (Я еще не решил).

Как я могу проверить свои результаты, полученные различными методами? Я пытался решить E(X) математически, но не смог найти какую-либо близкую форму. Эта функция может быть разделена на разные части, как

абсолютная (x)* двойная экспоненциальная плотность + другая функция, использующая более высокую степень (4) от X в обратной форме. Из-за абсолютного (x) и диапазона [-Inf, Inf] мы всегда должны делить его на [-Inf, 0] и [0, Inf]. Благодаря интеграции по частям я смог увидеть первую часть как (абсолютная (x) + (x^2/2) в бесконечном диапазоне) + Интеграл этой части не может быть найден математически.

Поэтому я использую следующий код, чтобы получить результат численного интегрирования как

Library(stats)
integrand <- function(x) {x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))}
integrate(integrand, lower = -Inf, upper = Inf)

таким образом, результат E(X)= 88,85766 с абсолютной ошибкой < 0,004

Результаты, которые я получаю с помощью этих методов, например, не похожи

(i) С помощью метода Симпсонов я получил E(X) = 0,3222642 и E(X^3)=677,0711.

simpson_v2 <- function(fun, a, b, n=100) {
    # numerical integral using Simpson's rule
    # assume a < b and n is an even positive integer
    if (a == -Inf & b == Inf) {
        f <- function(t) (fun((1-t)/t) + fun((t-1)/t))/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a == -Inf & b != Inf) {
        f <- function(t) fun(b-(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else if (a != -Inf & b == Inf) {
        f <- function(t) fun(a+(1-t)/t)/t^2
        s <- simpson_v2(f, 0, 1, n)
    } else {
        h <- (b-a)/n
        x <- seq(a, b, by=h)
        y <- fun(x)
        y[is.nan(y)]=0
        s <- y[1] + y[n+1] + 2*sum(y[seq(2,n,by=2)]) + 4 *sum(y[seq(3,n-1, by=2)])
        s <- s*h/3
    }
    return(s)
}

EX  <- function(x) x*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX, -Inf, Inf, n=100)

EX3 <- function(x) (x^3)*(abs(x)* exp(-0.5*abs(x))+(1/(1+(x-40)^4)))
simpson_v2(EX3, -Inf, Inf, n=100) 

(ii) Выборка по важности. Плотность моего предложения является нормальной со средним значением = 0 и стандартным отклонением =4. Резюме процесса выборки Важность, который я применяю, выглядит следующим образом

Предположим, я не могу сэмплировать из f (x), что верно, поскольку у него нет общеизвестной формы и нет встроенной функции в R, которую можно использовать для сэмплирования. Итак, я предлагаю другое распределение каркасно-канавных хвостов N(0, 4), чтобы брать выборки так, чтобы вместо оценки E(x) я оценивал E(x*f(x)/N(0,1)). Я использую следующий код для этого, который берет 100000 образцов из N (0,4)

X <- rnorm(1e5, sd=4)
Y <- (X)*(abs(X)*exp(-0.5*abs(X))+(1/(1+(x-40)^4)))/(dnorm(X, sd=4))
mean(Y)

Так как этот код требует случайной выборки из нормального распределения, поэтому каждый раз я получал разные ответы, но это что-то около -0.1710694, что почти похоже на 0.3222642. Я получил от метода Симпсонов. Но эти результаты очень отличаются от E(X)= 88,85766 от integrate (). Обратите внимание, что integrate () использует метод адаптивной квадратуры. Отличается ли этот метод от выборки Симпсонов и Важностей? Какого сходства в результатах я должен ожидать при сравнении этих методов?

Answer 1

Первый, EX а также EX3 определение неверно, вы пропустите минус под экспонентой

Ну вот несколько упрощений

• Если вы интегрируете эту часть x*(abs(x)* exp(-0.5*abs(x))), из -infinity... результат бесконечности будет 0

• Если вы интегрируете эту часть x^3*(abs(x)* exp(-0.5*abs(x))), из -infinity... результат бесконечности будет 0

• интеграл x^3/(1+(x-40)^4) из -infinity... бесконечность была бы бесконечностью, я рискну, вы получите бесконечные бесконечные логарифмы, см. http://integrals.wolfram.com/index.jsp?expr=%28xx x% 29% 2F% 281% + 2B+%28x-а% 29%29%5E4 & случайная = ложь

• интеграл x/(1+(x-40)^4) выглядит как нечто, похожее на обратную касательную, хотя онлайн-интегратор выдает уродливый вывод http://integrals.wolfram.com/index.jsp?expr=x%2F%281+%2B+%28x-a%29%5E4%29&random=false

ОБНОВИТЬ

Выглядит как твой EX было бы 40*\pi / \sqrt{2}

А также EX3 это не бесконечность, я могу ошибаться здесь

ОБНОВЛЕНИЕ 2

Ага, EX3 конечно, должно быть a^2*EX + \pi*a*3/\sqrt{2}, где a равно 40

ОБНОВЛЕНИЕ 3

Как уже отмечалось, для получения истинных значений требуется нормализация. EX а также EX3

N = 8 + \pi/\sqrt{2}

Вычисленные интегралы должны быть разделены на N чтобы получить подходящие моменты.