Search code examples
rnumerical-integration

Compute multiple Integral and plot them (with R)

I'm having trouble to compute and then plot multiple integral. It would be great if you could help me.

So I have this function

> f = function(x, mu = 30, s = 12){dnorm(x, mu, s)}

which i want to integrate multiple time between z(1:100) to +Inf to plot that with x=z and y = auc :

> auc = Integrate(f, z, Inf)

R return :

Warning message:
In if (is.finite(lower)) { :
  the condition has length > 1 and only the first element will be used

I have tested to do a loop :

while(z < 100){
z = 1
auc = integrate(f,z,Inf)
z = z+1}

Doesn't work either ... don't know what to do

Solution

  • Your problem is actually somewhat subtle, and in a certain sense gets to the core of how R works, so here is a slightly longer explanation.

    R is a "vectorized" language, which means that just about everything works on vectors. If I have 2 vectors A and B, then A+B is the element-by-element sum of A and B. Nearly all R functions work this way also. If X is a vector, then Y <- exp(X) is also a vector, where each element of Y is the exponential of the corresponding element of X.

    The function integrate(...) is one of the few functions in R that is not vectorized. So when you write:

    f   <- function(x, mu = 30, s = 12){dnorm(x, mu, s)}
auc <- integrate(f, z, Inf)

    the integrate(...) function does not know what to do with z when it is a vector. So it takes the first element and complains. Hence the warning message.

    There is a special function in R, Vectorize(...) that turns scalar functions into vectorized functions. You would use it this way:

    f   <- function(x, mu = 30, s = 12){dnorm(x, mu, s)}
auc <- Vectorize(function(z) integrate(f,z,Inf)$value)
z   <- 1:100
plot(z,auc(z), type="l")   # plot lines