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):
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
R:
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
...
The Gamma function of 0 or a negative integer is undefined, but it tends to +/-infinity as you approach those values. R returns NaN. It appears Python returns Inf, which works since the sign of the infinity doesn't appear to matter in your expression. To duplicate the Python behavioiur, you can replace gamma() with this function:
gamma <- function(x) if (round(x) == x && x <= 0) Inf else base::gamma(x)
Another approach would be to change your kum_mom() function so that it treats temp as zero when is.nan(temp) is TRUE, though this will give you warnings that you might want to suppress.