2016-01-23 44 views
0

我想将函数M2_11(如下所示)整合到x上,对于固定的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) 

TRUE

这是怎么回事检查结果?

+3

为什么'integrand = sapply(x,fM2_11)'而不是'fM2_11(x)'?我没有深入了解代码,但我想'fM2_11'可能不是矢量化的,而“集成”要求它是。阅读“集成”。 – nicola

+0

谢谢尼古拉!我欠你一个人情! –

回答

1

非常感谢尼古拉。问题已经解决了。

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

0.1588699与绝对误差< 1.7E-07

sum(integrand)*0.001 

0.1588705

不要期待回答如此简单!

+1

只是FYI,'Vectorize()'几乎就是'mapply()'的包装。在单一参数的情况下,单个结果函数将相当于'sapply()' –