Hi everyone,

I have a task to analytically compute a Cumulative distribution function (CDF) for the Gamma process with random effects. I developed the computation based on follows:

In my case, I set:

```
using Distributions, Statistics
x = 1:4000
α, β, δ, γ = 0.01252, 4.38e-3, 43.8, 0.0001
FT = 0.2775 - 0.1803
# compute cdf analytically
cdf_re_gp(t) = ccdf(FDist(2α*t, 2δ), FT/(δ*α*t*γ))
```

And I compare the above results with the Monte Carlo simulation, then I have the results:

As can be seen, there is some difference between these two results. For the Monte Carlo

simulation, I already take 11,000 trajectories. How can I make these two curves perfectly fit?

Thank you very much!

Just a guess, did you account for Distributions using the `Gamma(k, θ)`

instead of the `Gamma(α, β)`

parameterization?

No. I am from an engineering background and I did not check for other parameterizations (considering that they are equivalent).

Can you show a small demo for this? Thank you very much!

No demo necessary; it’s about what the arguments to `Gamma`

mean. The wikipedia page explains the 2 parameterizations: Gamma distribution - Wikipedia. Unfortunately papers and docstrings often don’t explain whether they use the shape-rate or shape-scale parameterization. The screenshot you shared uses \alpha and \beta for the parameters of the first Gamma distribution, which hints that they may be using shape-rate. But Distributions.jl docs (Univariate Distributions · Distributions.jl) indicate they use shape-scale.

Beyond that, there’s not much one can do to help without a complete code example.

1 Like

@sethaxen Thank you very much for the explanation.

I see that the two parametrizations only has a difference in the scale or the rate parameters.

Below is the Monte Carlo code I used to gather empirical failure times.

```
α, β, δ, γ = 0.01252, 4.38e-3, 43.8, 0.0001
function sim_hitting(r0, α0, β0, δ, γ, nbH; RE=false)
dt = 1
FT = 0.2775 - 0.1803
tins = []
rins = []
R_T1 = []
Trs::Vector{Float64} = []
Rtraj = []
for h in 1:nbH
r = r0
nbPM, nbCM, Ctot = 0.0, 0.0, 0.0
tsim = 0.0
rs = []
# add random effects
if RE
α = α0
β = rand(Gamma(δ, γ)) # 4.38e-3
else GP
α = α0
β = β0
end
while r < FT
# println(r)
# do PM at each inspection time
r += rand(Gamma(α*dt,β))
push!(rs, r)
tsim += dt
end
push!(Trs, tsim)
push!(Rtraj, [0; rs])
end
return Trs
end
```