2

Hi I am stuck with numerical integration of a function. I have this function:

Nd_f <- function(a,t) {
theta(t-a)*exp(-l*(1-exp(-a))) }  

which is used in another function defined below:

Nd <- function(s) {
integrate(Nd_f, lower = 0, upper = s, t=s)$value }  

where, theta() is a known function. So, using these functions I can evaluate Nd(t). But when I try to plot it using:

plot(Nd(0:500), log="y")

I get the following error:

Error in integrate(Nd_theta, lower = 0, upper = s, t = s) : evaluation of function gave a result of wrong length

I don't understand if I can evaluate it for all the values of t, why I cant plot it?

l=0.025, v= 0.001 and

theta <- function(t) { exp(-v*t) }

Thanks in advance!

0

1 Answer 1

6

I am assuming theta = exp and l = r = 1, so:

Nd_f <- function(a,t) exp(t-a)*exp(-(1-exp(-r*a)))

Function Nd <- function(s) integrate(Nd_f, lower = 0, upper = s, t=s)$value aims to evaluate the integral:

enter image description here

Note that integrate is not a vectorized function. A vectorized function can take vector input, and return a vector. For example, the exp function is vectorized, because:

exp(1:3)
# [1]  2.718282  7.389056 20.085537

But integrate is not. You are only allowed to pass in a scalar for lower and upper. So, there is no problem if you do:

Nd(1)
# [1] 1.273614

but it does not work when you do:

Nd(1:2)
# [1] 2.286086
# There were 15 or more warnings (use warnings() to see the first 15)
# warnings()
# Warning messages:
# 1: In t - a : longer object length is not a multiple of shorter object length

You need to wrap up your scalar function Nd to get a vectorized function. If you are really new to R, you might use a for loop:

Nd_vectorized_for <- function(s) {
  result <- numeric(length(s))
  for (i in 1:length(s)) {
    result[i] <- Nd(s[i])
    }
  result  ## or `return(result)`
  }

Now this function can take vector input and return a vector:

Nd_vectorized_for(1:2)
# [1] 1.273614 4.276839

People more experienced with R will suggest replacing for loop with *apply family function (read ?sapply to see this family):

Nd_vectorized_sapply <- function(s) sapply(s, Nd)

Nd_vectorized_sapply(1:2)
# [1] 1.273614 4.276839

But integrate is not a cheap operation, so there is no performance gain from sapply:

system.time(Nd_vectorized_for(sample(1:10,100000,replace=TRUE)))
#   user  system elapsed 
#  6.256   0.004   6.268 
system.time(Nd_vectorized_sapply(sample(1:10,100000,replace=TRUE)))
#   user  system elapsed 
#  6.200   0.004   6.212 

With a vectorized function, you can produce the plot you want:

plot(Nd_vectorized_for(1:50), log = "y")

enter image description here

3
  • 2
    Can't you use Vectorize i.e. Vectorize(Nd)(1:2) #[1] 1.273614 4.276839
    – akrun
    Jul 3, 2016 at 15:04
  • 1
    It is from experience.
    – akrun
    Jul 3, 2016 at 15:07
  • 1
    @akrun This was helpful, especially for multivariable arguments integration. Thanks!
    – VitalSigns
    Aug 26, 2016 at 13:00

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.