rfunctionsurvival-analysismaximize

Find the maximum of the function in R

I have the following function.

enter image description here

Let F(.) is the cumulative distribution function of the gamma distribution with shape = 1 and rate =1. The denominator is the survival function S(X) = 1 - F(X). The g(x) is the mean residual life function.

I wrote the following function in r.

x = 5
denominator = 1 -pgamma(x, 1, 1)
numerator = function(t) (1 - pgamma(t, 1, 1))

intnum  = integrate(numerator , x, Inf)

frac = intnum$value/denominator
frac

How can I find the maximum of the function g(x) for all possible values of X >= 0? Am I able to do this in r? Thank you very much for your help.

Solution

  • Before start, I defined the function you made

    surviveFunction<-function(x){
  denominator = 1 -pgamma(x, 1, 1)
  numerator = function(t) (1 - pgamma(t, 1, 1))

  # I used sapply to get even vector x
  intnum  = sapply(x,function(x){integrate(numerator , x, Inf)$value})
  
  frac = intnum/denominator
  return(frac)
}

    Then let's fit our function to function called 'curve' it will draw the plot with continuous data.

    The result is shown below:

    df = curve(surviveFunction, from=0, to=45)
plot(df, type='l')

    enter image description here And adjust the xlim to find the maximum value

    df = curve(surviveFunction, from=0, to=45,xlim = c(30,40))
plot(df, type='l')

    enter image description here

    And now we can guess the global maximum is located in near 35

    I suggest two options to find the global maximum.

    First using the df data to find maximum:

    > max(df$y,na.rm = TRUE)
 1.054248 #maximum value

> df$x[which(df$y==(max(df$y,na.rm = TRUE)))]
 35.55 #maximum value of x

    Second using the optimize:

    > optimize(surviveFunction, interval=c(34, 36), maximum=TRUE)

$maximum
[1] 35.48536

$objective
[1] 1.085282

    But the optimize function finds the not the global maximum value i think.

    If you see below

    optimize(surviveFunction, interval=c(0, 36), maximum=TRUE)

$maximum
[1] 11.11381

$objective
[1] 0.9999887

    Above result is not the global maximum I guess it is local maximum.

    So, I suggest you using first solution.