Используйте квадратуру Гаусса-Лагерра, чтобы приблизить интеграл в R - PullRequest
Новая
       76

Используйте квадратуру Гаусса-Лагерра, чтобы приблизить интеграл в R

0 голосов
/ 15 мая 2018

Я использую квадратуру Гаусса-Лагерра для аппроксимации следующего интеграла

First formula

Я написал функцию в R 

int.gl<-function(x)
{
 x^(-0.2)/(x+100)*exp(-100/x)
}

Сначала я использовалintegrate() функция для получения «истинного» значения, которое я дважды проверяю с помощью Wolfram Alpha 

integrate(int.gl, lower = 0, upper = Inf)$value
[1] 1.627777

Затем я использовал glaguerre.quadrature() function 

rule<-glaguerre.quadrature.rules(64, alpha = 0, normalized = F)[[64]]
glaguerre.quadrature(int.gl, lower = 0, upper = Inf, rule = rule, weighted = F)
[1] 0.03610346

Очевидно, что результатдалеко от истинного значения.В то время я думал, что должен был преобразовать эту функцию, чтобы включить термин ae ^ (- x).Итак, пусть X = 1 / Y.Я получил другую формулу, но тот же интеграл

Second formula

Таким же образом я использовал следующий код R 

int.gl2<-function(x)
{
 x^(-0.8)/(1+100*x)*exp(-100*x)
}
integrate(int.gl2, lower = 0, upper = Inf)$value
[1] 1.627777
glaguerre.quadrature(int.gl2, lower = 0, upper = Inf, rule = rule, weighted = F)
[1] 0.03937068

Ну, два разных значения.Есть ли у квадратуры Гаусса-Лагерра проблемы с нахождением этого интеграла?Есть ли какая-либо другая квадратура гауссовского типа, которая может помочь найти этот интеграл?

Примечание: Мне нужно использовать квадратурный тип Гаусса, так как я пытаюсь найти MLE некоторых параметров для настраиваемого распределения,Для простоты я просто зафиксировал эти параметры (эти константы в интеграле).Однако функция integrate() кажется менее «устойчивой», чем функция glaguerre.quadrature().(integrate() возвращает дивергентную ошибку при оптимизации логарифмической вероятности).

РЕДАКТИРОВАТЬ 1

Согласно Гансу.Комментарий W, я проверяю использование glaguerre.quadrature в следующем примере.Предположим, мы хотим найти следующий интеграл 

int.gl3<-function(x)
{
 ((x+100)/(x+200))^0.2*exp(-5*x)
}

glaguerre.quadrature(int.gl3, lower = 0, upper = Inf, rule = rule, weighted = F)
[1] 0.1741448
integrate(int.gl3, lower = 0, upper = Inf)$value
[1] 0.1741448

Кажется, использование glaguerre.quadrature является правильным.

Давайте теперь проверим преобразование.Я преобразую этот интеграл, задав 100Y = X, чтобы он включал в себя exp (-x). 

int.gl4<-function(x)
{
 0.01^0.2*x^(-0.8)/(1+x)*exp(-x)
}
integrate(int.gl4, lower = 0, upper = Inf)$value
[1] 1.627777
glaguerre.quadrature(int.gl4, lower = 0, upper = Inf, rule = rule, weighted = F)
[1] 0.9621667

Результат более близок к истинному значению, но не совсем равен.

РЕДАКТИРОВАТЬ 2

Вот мой полный пример

Интеграция и ложная сходимость оптимизации в R

...