Rинтегрировать: возвращает неправильное решение (использует неправильные квадратурные точки?)

У меня есть функция в R, которую я пытаюсь интегрировать, но для некоторых (экстремальных) значений параметров функции integrate возвращает неверное решение. Я полагаю, что проблема может заключаться в том, что integrate выбирает неправильные квадратурные точки для некоторых из этих экстремальных значений, но сначала я продемонстрирую проблему.

Функция, которую я хочу интегрировать, заключается в следующем.

integrandFunc_F <- function(x, func_u, func_u_lowerBar, 
  func_u_upperBar, func_mean_v, func_sigma_v, func_sigma_epsilon, 
  func_sigma_y, func_gamma, func_rho) {
#print(x);
p <- 1 - pnorm(func_u_upperBar,x,func_sigma_y);
q <- pnorm(func_u_lowerBar,x,func_sigma_y);
p <- p*(1-func_rho); q <- q*(1-func_rho);
alpha <- ifelse(func_gamma*(p+q) == 0, 0, pmax((func_gamma*p-q)/(func_gamma*(p+q)), 0));
g <- ifelse(x > func_u, dnorm(x,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))/(1-pnorm(func_u,func_mean_v,sqrt(func_sigma_v^2 + func_sigma_epsilon^2))), 0);
output <- alpha*g;
output
}

Когда я пытаюсь вычислить следующее, я получаю правильное решение 1:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 30, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

Однако, когда я пытаюсь вычислить следующее, я получаю неправильное решение 0:

integrate(integrandFunc_F, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
  func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
  func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)

Выше я указал, что, по моему мнению, проблема может быть связана с выбором точек квадратуры. Если вы раскомментируете #print(x) в приведенной выше функции, вы увидите, что в случае func_mean_v = 30 integrate устанавливает точки квадратуры, которые относительно велики/около 30. Однако в случае func_mean_v=50 после нескольких итераций integrate выбирает точки квадратуры, близкие к 0. Точки квадратуры около 0 не подходят для оценки этой функции, которая включает в себя нормальное распределение со средним значением func_mean_v..

Любые идеи о том, как решить эту проблему? Почему integrate в некоторых случаях выполняет итерацию по точкам квадратуры, близким к 0? Обратите внимание, что выбор func_mean_v = 30 и func_mean_v = 50, по общему признанию, является экстремальным параметром для этой функции, однако мне нужно уметь правильно рассчитывать такие случаи.


person Adam    schedule 30.08.2014    source источник
comment
Вы пробовали снизить допуск на сходимость интегрировать? Не существует общего правила для адаптивного численного интегрирования, которое будет работать для всех подынтегральных выражений (как узнать, где подынтегральное выражение не равно нулю между -Inf и +Inf?), поэтому иногда вы нужно помочь.   -  person baptiste    schedule 30.08.2014
comment
Я пробовал уменьшить rel.tol и увеличить количество subdivisions. Ни один из них не устранил проблему. Есть ли способ предоставить integrate набор точек квадратуры для начала? Для нормального распределения нетрудно узнать, где подынтегральная функция отлична от нуля.   -  person Adam    schedule 30.08.2014
comment
если вы знаете, где подынтегральная функция отлична от нуля, смещение и необязательное масштабирование переменной интегрирования - это то, что нужно имхо; это очень поможет integrate.   -  person baptiste    schedule 30.08.2014
comment
integrate() не может использовать ваши собственные квадратурные точки: он использует определенное квадратурное правило (Гаусса), поэтому веса подходят только для соответствующих узлов.   -  person baptiste    schedule 30.08.2014


Ответы (1)


вы можете сместить переменную интегрирования, чтобы центрировать пик,

wrapper <- function(x, func_mean_v, ...)
   integrandFunc_F(x+func_mean_v, func_mean_v=func_mean_v, ...)


integrate(wrapper, rel.tol = 1e-8, lower=-Inf, upper=Inf, func_u= 8, func_u_lowerBar= 8, 
          func_u_upperBar= 8, func_mean_v= 50, func_sigma_v= .1, func_sigma_epsilon= 2, 
          func_sigma_y= 1, func_gamma= 1/1.1, func_rho= .05)
# 1 with absolute error < 1.3e-09
person baptiste    schedule 30.08.2014
comment
Спасибо @baptiste. Можете ли вы отредактировать, чтобы показать точный код, который дал вам решение 1? Является ли оболочка функцией, внутри которой определен integrandFunc_F? Или они определяются отдельно? - person Adam; 30.08.2014
comment
это это точный код! вас может смутить отсутствие {} для оболочки, они не нужны для однострочников. - person baptiste; 30.08.2014
comment
Спасибо! Поэтому я сохраняю свое первоначальное определение integrandFunc_F, а затем дополнительно определяю функцию wrapper, как вы сделали здесь. - person Adam; 30.08.2014