approximated modified bessel second kind log

approximated_modified_bessel_second_kind_log.R

approximated_modified_bessel_second_kind_log = function(z, v, s = max(0,v-10)){

# s = 0 not approximated

  0.5 * ( log(pi) - log(2) ) + (-z - 0.5 * log(z)) +


  foreach(j = s:floor(v-0.5), .combine = c) %do% {

    lgamma(j + abs(v) - 0.5 + 1) -

    ( lgamma(j + 1) + lgamma(-j + abs(v) - 0.5 + 1) ) -

    j*(log(2) + log(z))

  } %>% matrixStats::logSumExp()
}