How to define a function value at a point for an integral in R?

Question

I would like to calculate the following double integral, with lower bound = -Inf and upper bound = Inf for both integrals. How can I define the function value to be zero for M=0 and/or leave out the numerical integration over M=0? It is a density function and I want to define the area over M=0 as zero when calculating probabilities.

The function can be seen here:

The R code is here, where one can see the "jump" due to my problem:

mum<-5.16
mub<-1.5
mur<-2.764
sm<-3.37
sb<-0.2
sr<-0.056

dk<-function(k){
integrate(function(r) { 
sapply(r, function(r) {
1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))},lower=Inf, upper=Inf)$value
})
}, lower=Inf, upper=Inf)$value

}

curve(sapply(x,dk), from=-2, to=18, lwd=2,col="red")

your code does not run as is. There is an unexpected '}'. Also, no data provided to test the function for your desired behaviour. Could you edit the question to address these issues. I expect that replacing your integrate function with function(m) { if (m==0) 0 else 1/abs(m)*exp(-(m-m3)^2-((k-r)/m-m2)^2) } will address your problem, but need data to be sure — dww, May 09 '16 at 14:36
PLease see http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example for how to create a MWE. — dww, May 09 '16 at 14:37
Thanks for your input, I corrected the brackets. However, I have no data to fit, it is just parameters which I have to set. The function I use in real and the code is much longer. For presenting the problem of dividing throug m I simplified it. It is a continous distribution function and I can define the area over point as zero, which I need to implement in R. — Seb, May 09 '16 at 14:45
I cannot use if(m==0) 0 else because m is not an object, it is the integration variable. — Seb, May 09 '16 at 14:54

score 2 · Answer 1 · answered May 09 '16 at 16:13

One approach would be to simply not to evaluate your function dk at values where its value is invalid as in e.g.

mum<-5.16
mub<-1.5
mur<-2.764
sm<-3.37
sb<-0.2
sr<-0.056

  dk<-function(k){
    integrate(function(r) { 
      sapply(r, function(r) {
        1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))}, lower=Inf, upper=Inf)$value
      })
    }, lower=Inf, upper=Inf)$value
  }

x <- c(seq(-2,-0.1, by=0.2), seq(0.1,18,by=0.2))
y <- sapply(x,dk)
plot(y~x, type="l", col="red")

The other way is to have your function dk return NA at values for which the function fails as in e.g.

  dk<-function(k){
    if ((k-mur)<0.4 & (k-mur)>-0.4) NA
    else {
    integrate(function(r) { 
      sapply(r, function(r) {
        1/(2*pi*sr^2)^0.5*exp(-(r-mur)^2/(2*sr^2))*integrate(function(m) {1/(abs(m)*2*pi*sm*sb)*exp(-(m-mum)^2/(2*sm^2)-((k-r)/m-mub)^2/(2*sb^2))}, lower=Inf, upper=Inf)$value
      })
    }, lower=Inf, upper=Inf)$value
    }
  }

curve(sapply(x,dk), from=-2, to=18, lwd=2,col="red")

How to define a function value at a point for an integral in R?

1 Answers1