integrate()에 점프합니까?

Question

몇 가지 적분을 계산해야합니다. 나는 R이 이것을하기위한 올바른 소프트웨어가 아니라는 것을 알고있다. 그러나 R에서 모든 것을 다했듯이 적분은 단지 1 차원 적이기 때문에 문제가 없을 것이라고 생각했다.

어쨌든 통계 패키지의 integrate() 함수를 사용하면 적분의 상한이 증가 할 때 점프가 발생합니다. 여기 플롯이다 : 나는, 적분 대신 합을 구축하거나 cubature 패키지에서 기능 adaptIntegrate()를 사용하는 경우에

은, 그러나, 결과는 다음과 같습니다

이 코드를 복제하려면 다음 코드를 작성하십시오. 아마 더 쉬운 예가 있다는 것을 압니다. 그러나 이것은 실제로 내가 직면하고있는 경우의 1 : 1입니다. integrate() 사용 :

v_p = 11269 
d_p = seq(0, v_p, 100) 
tau = 0.35 
interest = 0.05 
time_t = 1 
mu = 0 
sigma_p = 0.1099313 
nu_p = 0.0101552 
ts = c() 
bc = c() 
for (i in 1:length(d_p)){ 
    ts = c(ts,integrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$value) 
    bc = c(bc,integrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$value) 
} 
shv = ts + bc 
plot(d_p,shv, main = 290801) 
abline(v = 9079, col="red") 

정확히 동일한 코드를, 단지 adaptIntegrate()를 사용하여 :

는

v_p = 11269 
d_p = seq(0, v_p, 100) 
tau = 0.35 
interest = 0.05 
time_t = 1 
mu = 0 
sigma_p = 0.1099313 
nu_p = 0.0101552 
ts = c() 
bc = c() 
for (i in 1:length(d_p)){ 
    ts = c(ts,adaptIntegrate(function(x) tau*(x-v_p-interest*d_p[i])*dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), v_p+interest*d_p[i], 10000000)$integral) 
    bc = c(bc,adaptIntegrate(function(y) nu_p * y * dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)), 0, d_p[i])$integral) 
} 
shv = ts + bc 
plot(d_p,shv, main = 290801) 
abline(v = 9079, col="red") 

는 아무도 왜 이런 일 아이디어가 있습니까?

Answer 1

당신은 ?integrate에 주의 두 번째 단락을 읽어야 할 사람 :

무한 간격에 걸쳐 통합 때문에 명시 적으로하지 않고 단지 끝점으로 많은 수를 사용하는 것보다. 이것은 정답의 기회를 증가시킵니다. 무한 구간에서의 적분이 유한 한 함수는 그 구간의 대부분 동안 0에 가까워 야합니다. 내가 Inf를에 10000000 당신의 상한을 변경하는 경우

(뿐만 아니라 읽을 가치가 주의 나머지.)

나는 합리적인 보이는 결과를 얻을. (그러나 이런 종류의 문제는 일 수 있습니다.은 일반적으로 완전히 해결됩니다 ... 특정 문제에 대해 간격, 시작점 수 등을 조정해야하는 경우가 있습니다.)

코드 코어 INF는 명확성을 위해 약간, 다시 대체 :

f1 <- function(x) tau*(x-v_p-interest*d_p[i])* 
     dlnorm(x, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)) 
f2 <- function(y) nu_p * y * 
     dlnorm(y, log(v_p)+(mu - 0.5 * sigma_p^2) * time_t, sigma_p*sqrt(time_t)) 
for (i in 1:length(d_p)){ 
    ts = c(ts,integrate(f1, v_p+interest*d_p[i], Inf)$value) 
    bc = c(bc,integrate(f2, 0, d_p[i])$value) 
} 

(그건 그렇고, 당신은 또한 R에서 개체를 성장하지 않도록해야한다 전체 벡터를 할당 첫번째보다는 반복적으로 그들에게 추가.)