Limited Expected Value for a Gamma Distribution

A student recently asked me how to write a script for this function, that involves the incomplete gamma function. Here is the R version

# E[X^x] for the gamma function

library(tidyverse)

elv <- function(shape, scale, lim) {
  
  gamma_pdf <- function(x, shape, scale) dgamma(x, shape = shape, scale = scale)  
  int <- integrate(function(x) x * gamma_pdf(x, shape, scale), 0, lim)$value  
  cdf_lim <- pgamma(lim, shape = shape, scale = scale)  
  return(int + lim * (1 - cdf_lim))  
}

# Example
shape = 2; scale = 10
df <- data.frame(
  deductible = seq(0, 10, 0.1),
  ELV = sapply(seq(0, 10, 0.1), FUN = elv, shape = shape, scale = scale)
)