So, I’m trying to compute the mean of a distribution, it goes smoothly in python and the results match the article I’m reading, but it doesn’t work in R, probably because of the negative inputs for the gamma function since b-i is going to be negative most of the time. Is there any package or modification I could do to make it works? The resources that I have for this distribution are mostly in R, so I don’t really need it in python, just saw it was working there.

Python:

import math
import numpy as np
from scipy.special import gamma
def kum_mom(lambd, b, beta, r):
  limite = 100
  i = 0
  soma = 0
  while i <= limite:
    temp=(((-1)**i)*lambd*b*beta*gamma(b)*gamma(1-(r/beta)))/(math.factorial(i)*beta*((lambd*(i+1))**(1-              (r/beta)))*gamma(b-i))
    soma=soma+temp
    i = i + 1
  return soma

kum_mom(1,1,3,1)
1.3541179394264005

R:

library(pracma)
kum_mom <- function(lambda, b, beta, r) {
  limite <- 100
  i <- 0
  soma <- 0
  while (i <= limite) {
    temp <- ((-1)^i * lambda * b * beta * gamma(b) * gamma(1 - (r / beta))) / (factorial(i) * beta * ((lambda * (i + 1))^(1 - (r / beta))) * gamma(b - i))
    soma <- soma + temp
    i <- i + 1
  }
  return(soma)
}

kum_mom(1, 1, 3, 1)
[1] NaN
There were 50 or more warnings (use warnings() to see the first 50)

warnings()
1: In gamma(b - i) : NaNs produced
2: In gamma(b - i) : NaNs produced
3: In gamma(b - i) : NaNs produced
...