stable RNG

Author: PharmCat

.github/workflows/Documenter.yml

@@ -24,7 +24,7 @@ jobs:
       - uses: actions/checkout@v2
       - uses: julia-actions/setup-julia@latest
         with:
-          version: '1.6'
+          version: '1.5'
       - name: Install dependencies
         run: julia --project=docs/ -e 'using Pkg; Pkg.develop(PackageSpec(path=pwd())); Pkg.instantiate()'
       - name: Build and deploy

src/sim.jl

@@ -5,17 +5,17 @@
 Bioequivalence power simulation.
 
 """
-function besim(t::CTask{P, D, Bioequivalence, Power}; nsim::Int = 100, rsabe::Bool = false, rsabeconst = 0.760, seed = 0)  where P where D
-    if seed != 0  Random.seed!(seed) end
+function besim(t::CTask{P, D, Bioequivalence, Power}; nsim::Int = 100, rsabe::Bool = false, rsabeconst = 0.760, seed = 0, rng = MersenneTwister())  where P where D
+    if seed != 0  Random.seed!(rng, seed) end
     df      = t.design.df(t.objective.val)
     sef     = sediv(t.design, t.objective.val)
     CHSQ    = Chisq(df)
     tval    = quantile(TDist(df), 1 - t.hyp.alpha)
     σ̵ₓ      = sef * t.param.a.sd
     pow     = 0
     for i = 1:nsim
-        smean   = rand(ZDIST) * σ̵ₓ + (t.param.a.m - t.param.b.m)
-        σ²      = rand(CHSQ)  * t.param.a.sd  ^ 2 / df
+        smean   = rand(rng, ZDIST) * σ̵ₓ + (t.param.a.m - t.param.b.m)
+        σ²      = rand(rng, CHSQ)  * t.param.a.sd  ^ 2 / df
         σ       = sqrt(σ²)
         hw      = tval * σ * sef
         θ₁      = smean - hw
@@ -41,16 +41,16 @@ end
 Proportion difference power simulation.
 
 """
-function ctsim(t::CTask{DiffProportion{P1, P2}, D, H, Power}; nsim = 100, seed=0, method = :default, dropout = 0.0)  where P1 <: Proportion where P2 <: Proportion where D where H <: AbstractHypothesis
-    if seed != 0  Random.seed!(seed) end
+function ctsim(t::CTask{DiffProportion{P1, P2}, D, H, Power}; nsim = 100, method = :default, dropout = 0.0, seed = 0, rng = MersenneTwister())  where P1 <: Proportion where P2 <: Proportion where D where H <: AbstractHypothesis
+    if seed != 0  Random.seed!(rng, seed) end
     n₁      = getval(t.objective)
     n₂      = Int(ceil(getval(t.objective) / t.k))
     pow     = 0
     D₁      = Binomial(n₁, getval(t.param.a))
     D₂      = Binomial(n₂, getval(t.param.b))
     for i = 1:nsim
-        xn₁      = rand(D₁)
-        xn₂      = rand(D₂)
+        xn₁      = rand(rng, D₁)
+        xn₂      = rand(rng, D₂)
         if  checkhyp(t.hyp, DiffProportion(Proportion(xn₁, n₁), Proportion(xn₂, n₂)); method = method) pow += 1 end
     end
     return pow/nsim
@@ -122,15 +122,15 @@ end
 Means power simulation.
 
 """
-function ctsim(t::CTask{P, D, H, Power}; nsim = 100, seed = 0, method = :default, dropout = 0.0)  where P <: DiffMean where D where H
-    if seed != 0  Random.seed!(seed) end
+function ctsim(t::CTask{P, D, H, Power}; nsim = 100, method = :default, dropout = 0.0, seed = 0, rng = MersenneTwister())  where P <: DiffMean where D where H
+    if seed != 0  Random.seed!(rng, seed) end
     n₁      = getval(t.objective)
     n₂      = Int(ceil(getval(t.objective) / t.k))
     pow     = 0
 
     for i = 1:nsim
-        n1      = rand(ZDIST, n₁) .* t.param.a.sd .+ getval(t.param.a)
-        n2      = rand(ZDIST, n₂) .* t.param.b.sd .+ getval(t.param.b)
+        n1      = rand(rng, ZDIST, n₁) .* t.param.a.sd .+ getval(t.param.a)
+        n2      = rand(rng, ZDIST, n₂) .* t.param.b.sd .+ getval(t.param.b)
         if  checkhyp(t.hyp, DiffMean(Mean(n1), Mean(n2)); method = method) pow += 1 end
     end
     return pow/nsim

test/citest.jl

@@ -169,6 +169,7 @@ println(" ---------------------------------- ")
     @test ci.upper    ≈ Inf
     @test ci.estimate ≈ Inf
 
+
     #-- cli
 
     ci = ClinicalTrialUtilities.rrpropci(30, 100, 40, 90; alpha=0.05, method=:cli)

test/sim.jl

@@ -2,13 +2,13 @@ println(" ---------------------------------- ")
 @testset "#14 Simulations         " begin
 
 t      = ClinicalTrialUtilities.bepower(cv=0.2, n=20).task
-result = ClinicalTrialUtilities.besim(t; nsim = 100, seed = 1234)
-@test result == 0.83
+result = ClinicalTrialUtilities.besim(t; nsim = 100, seed = 1234, rng = StableRNG(0))
+@test result == 0.82
 Base.show(io, t)
 
 t      = ClinicalTrialUtilities.ctask(param=:mean, hyp=:ns, group=:two, alpha=0.05, n = 108, sd=10, diff=5, a=7, b=0, k=1)
-result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234)
-@test result == 0.431
+result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234, rng = StableRNG(0))
+@test result == 0.439
 Base.show(io, t)
 #=
 julia> ClinicalTrialUtilities.ctpower(t)
@@ -28,8 +28,8 @@ Power: 0.431398
 =#
 
  t      = ClinicalTrialUtilities.ctask(param=:mean, hyp=:ea, group=:two, alpha=0.05, n = 108, sd=5, diff=0, a=5, b=6, k=1)
- result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234)
- @test result == 0.315
+ result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234, rng = StableRNG(0))
+ @test result == 0.323
  Base.show(io, t)
 
 #=
@@ -50,8 +50,8 @@ Power: 0.312274
 =#
 
 t      = ClinicalTrialUtilities.ctask(param=:mean, hyp=:ei, group=:two, alpha=0.05, n = 108, sd=10, diff=5, a=5, b=4, k=1)
-result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234)
-@test result == 0.885
+result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234, rng = StableRNG(0))
+@test result == 0.884
 Base.show(io, t)
 #=
 julia> ClinicalTrialUtilities.ctpower(t)
@@ -77,8 +77,8 @@ ClinicalTrialUtilities.DiffProportion(ClinicalTrialUtilities.Proportion(30, 100)
 ClinicalTrialUtilities.Parallel(),
 ClinicalTrialUtilities.Superiority(-0.15, -0.15, 0.05),
 ClinicalTrialUtilities.Power(100), 1.0)
-result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, method = :nhs, seed = 1234)
-@test result == 0.169
+result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, method = :nhs, seed = 1234, rng = StableRNG(0))
+@test result == 0.214
 Base.show(io, t)
 
 #=
@@ -103,8 +103,8 @@ ClinicalTrialUtilities.DiffProportion(ClinicalTrialUtilities.Proportion(75, 100)
 ClinicalTrialUtilities.Parallel(),
 ClinicalTrialUtilities.Equivalence(-0.2, 0.2, 0.05),
 ClinicalTrialUtilities.Power(132), 1.0)
-result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234)
-@test result == 0.894
+result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, seed = 1234, rng = StableRNG(0))
+@test result == 0.897
 Base.show(io, t)
 
 #=
@@ -129,8 +129,8 @@ ClinicalTrialUtilities.DiffProportion(ClinicalTrialUtilities.Proportion(30, 100)
 ClinicalTrialUtilities.Parallel(),
 ClinicalTrialUtilities.Equality(0, 0.05),
 ClinicalTrialUtilities.Power(100), 1.0)
-result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, method = :mn, seed = 1234)
-@test result == 0.343
+result = ClinicalTrialUtilities.ctsim(t; nsim = 1000, method = :mn, seed = 1234, rng = StableRNG(0))
+@test result == 0.315
 Base.show(io, t)
 
 #=