Борьба с гамма-функцией

Итак, я пытаюсь вычислить среднее значение распределения, в Python оно проходит гладко, и результаты соответствуют статье, которую я читаю, но это не работает в R, вероятно, из-за отрицательных входных данных для гамма-функции, поскольку b-i большую часть времени будет отрицательным. Есть ли какой-нибудь пакет или модификация, которую я мог бы сделать, чтобы он работал? Ресурсы, которые у меня есть для этого дистрибутива, в основном написаны на R, поэтому на Python они мне особо не нужны, просто увидел, что там они работают.

Питон:

import math
import numpy as np
from scipy.special import gamma
def kum_mom(lambd, b, beta, r):
  limit = 100
  i = 0
  sum = 0
  while i <= limit:
    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))
    sum = sum + temp
    i = i + 1
  return sum

kum_mom(1, 1, 3, 1)
1.3541179394264005

Р:

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

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
...