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

**URL:** <https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747>\
**Category:** Statistics\
**Tags:** distributions, cdf\
**Created:** [February 16, 2023, 9:50pm UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747 "2023-02-16T21:50:41Z")\
**Posts on this page:** 5\
**Page:** 1

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 16, 2023, 9:50pm UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747/1 "2023-02-16T21:50:41Z")

</div>

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:  
 ![re-gp](https://global.discourse-cdn.com/julialang/original/3X/c/6/c66023514be47074b6bf5b9712bd41e19be8446f.png)  
In my case, I set:

```julia
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:  
 ![re-gp1](https://global.discourse-cdn.com/julialang/original/3X/4/6/460a60174ce3a37207567a1cf167ae90244c31bb.png)

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!

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [February 20, 2023, 7:03am UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747/2 "2023-02-20T07:03:41Z")

</div>

Just a guess, did you account for Distributions using the `Gamma(k, θ)` instead of the `Gamma(α, β)` parameterization?

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [February 20, 2023, 7:13am UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747/3 "2023-02-20T07:13:39Z")

</div>

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!

---

<div class="post-metadata">

**Author:** ![sethaxen](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/sethaxen/32/35604_2.png) [@sethaxen](https://discourse.julialang.org/u/sethaxen)\
**Post date:** [March 8, 2023, 10:21am UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747/4 "2023-03-08T10:21:19Z")

</div>

No demo necessary; it’s about what the arguments to `Gamma` mean. The wikipedia page explains the 2 parameterizations: [Gamma distribution - Wikipedia](https://en.wikipedia.org/wiki/Gamma_distribution#Definitions). 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](https://juliastats.org/Distributions.jl/stable/univariate/#Distributions.Gamma)) indicate they use shape-scale.

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

---

<div class="post-metadata">

**Author:** ![Jian\_ZUO](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jian_zuo/32/33738_2.png) [@Jian\_ZUO](https://discourse.julialang.org/u/Jian_ZUO)\
**Post date:** [March 8, 2023, 4:36pm UTC](https://discourse.julialang.org/t/cumulative-distribution-function-results-are-not-perfectly-fit-for-simulation-and-analytical-methods/94747/5 "2023-03-08T16:36:47Z")

</div>

@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.

```julia
α, β, δ, γ = 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

```