Search code examples
rnumerical-integrationuniroot

Solve improper double integral using integrate and uniroot functions


We have a function. t ~ Weibull(alpha, lambda) and c ~ Exponential(beta):

enter image description here

Given p = 0.10, alpha = 1, lambda = 4. Find the value of beta.

We want to integrate this function for t then to c. Then find the value of beta where integral equals to p using uniroot function.

See the code below:

alpha = 1
lambda = 4
p = 0.10
func1 <- function(t, c, beta) {alpha * lambda * exp(-lambda * t^ alpha)*
                           beta * exp(- beta * c) }
func2 <- function(c, beta){integrate(func1, lower = c, upper = Inf, c=c, 
beta=beta)}
func3 <- function(beta){integrate(func2, lower = 0, upper = Inf, beta = 
beta)$value - cen.p}
uniroot(func3 ,lower = 0.001, upper = 10, extendInt = "yes")$root

However it throws the error:

Error in integrate(func1, lower = c, upper = Inf, c = c, beta = beta) : length(lower) == 1 not TRUE

Answer should be 0.444


Solution

  • I corrected typos (substituted cen.p to p) and vectorized function arguments for func2 and func3, since the integrate function returns one value (scalar). However as a first argument integrate should accept vector of numeric values, not a scalar.

    alpha <- 1
    lambda <- 4
    p <- 0.10
    func1 <- function(t, c, beta)
      alpha * lambda * t^(alpha - 1) * exp(-lambda * t^alpha) * beta * exp(-beta * c)
    
    func2 <- function(c, beta)
      integrate(func1, lower = c, upper = Inf, c = c, beta = beta)$value)
    
    func3 <- function(beta)
      integrate(Vectorize(func2), lower = 0, upper = Inf, beta = beta)$value - p
    
    uniroot(Vectorize(func3), lower = 0.001, upper = 10, extendInt = "yes")$root
    

    Output:

    [1] 0.4444242.