2016-01-23 3 views
0

Я хотел бы интегрировать функцию M2_11 (далее) по х при фиксированном theta = c(2,0.8), c = 1.1, a=c(1,1) и A = matrix(c(1/0.8,0.03,0.03,2/0.8),nrow=2,ncol=2).«интегрировать» функции в R дает ошибочные результаты

M2_11 = function(x, theta, c, a, A){ 
return((score1(x,theta)-a[1])^2* (weight(x, theta, c, a, A))^2 * f(x, theta)) 
} 

интегрирующий функция R дает следующие результаты

theta = c(2,0.8) 
c = 1.1 
a=c(1,1) 
A = matrix(c(1/0.8,0.03,0.03,2/0.8),nrow=2,ncol=2) 
integrate(M2_11, lower = 1e-100, upper = 10 ,subdivisions = 10000, theta,c,a,A) 

0,0006459957 с абсолютной погрешностью < 4.5e-05

Выполнение интеграции другой способ дает тот же результат

fM2_11 = function(x){M2_11(x,theta,c,a,A)} 
integrate(fM2_11, lower = 1e-100, upper = 10,subdivisions = 10000) 

0,0006459957 с абсолютной погрешностью < 4.5e-05

Результат, что функция интеграции дает, однако, явно не так:

x = seq(1e-100,10,by=0.001) 
integrand = sapply(x,fM2_11) 

enter image description here

площадь под кривая явно превышает 0,00066

I Также проверьте результат с помощью цикла

loop_result = rep(NA,length(x)) 
for (i in 1:length(x)){ 
    loop_result[i] = M2_11(x[i],theta,c,a,A) 
} 
table(integrand==loop_result) 

ИСТИНА

Что происходит?

+3

Почему 'подынтеграл = sapply (x, fM2_11)' вместо просто 'fM2_11 (x)'? Я не копал код, но я думаю, что 'fM2_11' может быть не векторизованным, а' integrate' требует его. Читайте '? Integrate'. – nicola

+0

Спасибо Никола! Я твой должник! –

ответ

1

спасибо Nicola. Проблема решена.

integrate(Vectorize(fM2_11), lower = 1e-100, upper = 10 ,subdivisions = 10000) 

0,1588699 с абсолютной погрешностью < 1.7E-07

sum(integrand)*0.001 

0,1588705

Никогда не ожидать ответа быть так просто!

+1

Just FYI, 'Vectorize()' в значительной степени просто оболочка для 'mapply()'. В случае единственного аргумента единственная функция результата будет эквивалентна 'sapply()' –

Смежные вопросы