2

I'm calculating the probability mass function for a count variable and the normalization term is an infinite sum of the form ∑f(n), where the sum goes over all non-negative integers (0-infinity). I'm looking for a function in R that approximates such sum. After some research I've found that classical procedures are approximations like the Laplace method for sums or the Euler-Maclaurin sum formula, but I'm unable to find a function for that in R. The function f(n) becomes decreasing after some n and converges to 0. The function is

f(n) = exp(an + b*sqrt(n) - ln(n!)),

where a and b are some constants.

gregorp
  • 255
  • 1
  • 3
  • 11

1 Answers1

2

Finding the infinite sum of values of a function is the same as integrating this function over the same interval. You can use function integrate():

#Define the function
f <- function(x){ 1/(x**2)}

# First check on a segment
lower = 0.5
upper = 2
sum(f(seq(lower, upper, by=0.0001)))
# [1] 15002.13

# Integrate
integrate(f, lower, upper)
# 1.5 with absolute error < 3.8e-09

# For upper boundary to be infinity:
lower = 0.5
upper = Inf
integrate(f, lower, upper)
# 2 with absolute error < 8.4e-11

If integration is over the integer values, you can try to approximate the upper boundary and calculate (though for some functions it might be very slow):

sum(f(1:(2^25)))
# [1] 1.644934

sum(f(1:(2^26))) # Check how much the value changes for even longer vector
# [1] 1.644934

Added after the function definition was added to the question: A few notes about this specific function:

f(n) = exp(an + b*sqrt(n) - ln(n!))

There is a term ln(n!) in this function that R will have a very hard time to estimate, since it will need to calculate N! for a very large number. So the best approach for this particular case is to find upper and lower boundary functions and use them for approximation.

Katia
  • 3,784
  • 1
  • 14
  • 27
  • But this only applies for very small fragments (which becomes the Riemann integral if the size of the fragments goes to 0). In my case the sum is over non-negative integers. – gregorp May 07 '18 at 13:56
  • I see I was not clear enough and wrote non-negative numbers instead of integers, I apologize. – gregorp May 07 '18 at 14:02
  • So what you're proposing is that I truncate the sum? As in sum only to some very high n? Just truncating would not be ok, because this is a normalization term and if I truncate, then the probabilities will sum to more than 1. Maybe I could find a continuous function as an upper bound above some n and then approximate that. – gregorp May 07 '18 at 15:57
  • 1
    Yes, my suggestion with trancation of the sum will only approximate the sum. I do not know any packages that would perform the sum over the integers, so the approach of finding an upper and lower limits for your functions will be the best approach... – Katia May 07 '18 at 16:06
  • Also , you might want to correct your function definition - you use"n" in f(n), but then use x within the function. – Katia May 07 '18 at 16:07