# Cumulative distribution function results are not perfectly fit for simulation and analytical methods

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 = []
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