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