Renewal Function for Weibull Distribution

The renewal function for Weibull distribution m(t) with t = 10 is given as below.

I want to find the value of m(t). I wrote the following r code to compute m(t)

last_term = NULL
gamma_k = NULL
n = 50
for(k in 1:n){
  gamma_k[k] = gamma(2*k + 1)/factorial(k)
}

for(j in 1: (n-1)){
  prev = gamma_k[n-j]
  last_term[j] = gamma(2*j + 1)/factorial(j)*prev
}

final_term = NULL
find_value = function(n){
  for(i in 2:n){
  final_term[i] = gamma_k[i] - sum(last_term[1:(i-1)])
  }
  return(final_term)
}
all_k = find_value(n)

af_sum = NULL
m_t = function(t){
for(k in 1:n){
af_sum[k] = (-1)^(k-1) * all_k[k] * t^(2*k)/gamma(2*k + 1)
}
  return(sum(na.omit(af_sum)))
}
m_t(20)

The output is m(t) = 2.670408e+93. Does my iteratvie procedure correct? Thanks.

I don't think it will work. First, lets move Γ(2k+1) from denominator of m(t) into A_k. Thus, A_k will behave roughly as 1/k!.

In the nominator of the m(t) terms there is t^2k, so roughly speaking you're computing sum with terms

100^k/k!

From Stirling formula

k! ~ k^k, making terms

(100/k)^k

so yes, they will start to decrease and converge to something but after 100th term

Anyway, here is the code, you could try to improve it, but it breaks at k~70

N <- 20
A <- rep(0, N)

# compute A_k/gamma(2k+1) terms
ps <- 0.0 # previous sum
A[1] = 1.0
for(k in 2:N) {
    ps <- ps + A[k-1]*gamma(2*(k-1) + 1)/factorial(k-1)
    A[k] <- 1.0/factorial(k) - ps/gamma(2*k+1)
}

print(A)

t <- 10.0
t2 <- t*t

r <- 0.0
for(k in 1:N){
    r <- r + (-t2)^k*A[k]
}

print(-r)

UPDATE

Ok, I calculated A_k as in your question, got the same answer. I want to estimate terms A_k/Γ(2k+1) from m(t), I believe it will be pretty much dominated by 1/k! term. To do that I made another array k!*A_k/Γ(2k+1), and it should be close to one.

Code

N <- 20
A <- rep(0.0, N)

psum <- function( pA, k ) {
    ps <- 0.0
    if (k >= 2) {
        jmax <- k - 1
        for(j in 1:jmax) {
            ps <- ps + (gamma(2*j+1)/factorial(j))*pA[k-j]
        }
    }
    ps
}

# compute A_k/gamma(2k+1) terms
A[1] = gamma(3)
for(k in 2:N) {
    A[k] <- gamma(2*k+1)/factorial(k) - psum(A, k)
}

print(A)

B <- rep(0.0, N)
for(k in 1:N) {
    B[k] <- (A[k]/gamma(2*k+1))*factorial(k)
}

print(B)

shows that

1. I got the same A_k values as you did.
2. B_k is indeed very close to 1

It means that term A_k/Γ(2k+1) could be replaced by 1/k! to get quick estimate of what we might get (with replacement)

m(t) ~= - Sum(k=1, k=Infinity) (-1)^k (t²)^k / k! = 1 - Sum(k=0, k=Infinity) (-t²)^k / k!

This is actually well-known sum and it is equal to exp() with negative argument (well, you have to add term for k=0)

m(t) ~= 1 - exp(-t²)

Conclusions

1. Approximate value is positive. Probably will stay positive after all, A_k/Γ(2k+1) is a bit different from 1/k!.

2. We're talking about 1 - exp(-100), which is 1-3.72*10^-44! And we're trying to compute it precisely summing and subtracting values on the order of 10¹⁰⁰ or even higher. Even with MPFR I don't think this is possible.

Another approach is needed