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.