@storopoli Much appreciated.
Is there a way to output the trace of parameters during the rapid exploration phase in SAEM? I want to look at the progression over the iterations to better understand the algorithm. Thank you.
Bob.
You can set verbosity in your Pumas.SAEM.
You will find all the options in the documentation from the help which is usually well populated for most Pumas functions from the REPL.
One can get the docs of any function in Julia using ? function. See below for SAEM
help?> Pumas.SAEM
SAEM(; iters::Tuple{Int,Int,Int} = (200,100,100), α::Float64 = 0.7, λ::Float64 = 0.975, repeat::Int = 1, verbosity::Int = 0)
• iters: A tuple indicating the number of iterations to run in the exploratory, middle, and smoothing phases of SAEM.
• α: The learning rate for the middle phase. On the kth update, α^k weight is placed on the kth iteration, and 1 - α^k on the previous.
• λ: Maximum standard deviation change rate. A value of 0.95 means that ω_new = clamp(ω_new, ω_old * 0.95, ω_old / 0.95) on each iteration. Slowing the
decliine in ωs and σs increases exploration, while limiting the growth prevents excessive flattening of the log likelihood.
• repeat: Indicates how many times to repeat iters. Defaults to 1.
• verbosity: How verbose should it be? 0: silent, 1: print on each repeat, 2: print once per iteration.
This verbosity printing may not be really useful without a visual/graphical representation over time. This is a feature in the pipeline and not available now.
@vijay Thank you. I was using the Documentation. repeat argument does not seem to be supported? The only line I added was the repeat=1 to check.
rngv = [MersenneTwister(1941964947i + 1) for i ∈ 1:Threads.nthreads()]
ka1cmt_em_fit = fit(ka1cmt_saem,
pop_nmle,
param_saem,
#Pumas.SAEM(),
Pumas.SAEM(iters=(5000,1000,10000)),
repeat = 1,
ensemblealg = EnsembleThreads(),
rng=rngv)
ERROR: LoadError: MethodError: no method matching fit(::Pumas.PumasEMModel{(1, 1, 1), 0, (), (), (), (:Depot, :Central), (:cp,), Nothing, Tuple{Pumas.Formula{:Ka, (1,), :LogNormal}, Pumas.Formula{:CL, (1,), :LogNormal}, Pumas.Formula{:Vc, (1,), :LogNormal}}, Nothing, var"#55#57", Depots1Central1, var"#56#58", Tuple{Pumas.Error{:CONC, :cp, :Normal}}, Nothing}, ::Vector{Subject{NamedTuple{(:CONC,), Tuple{Vector{Union{Missing, Float64}}}}, Pumas.ConstantCovar{NamedTuple{(:WEIGHTB, :SEX, :DOSE), Tuple{Float64, InlineStrings.String7, Int64}}}, Vector{Pumas.Event{Float64, Float64, Float64, Float64, Float64, Float64, Int64}}, Vector{Float64}}}, ::NamedTuple{(:Ka, :CL, :Vc), Tuple{Int64, Int64, Int64}}, ::Pumas.SAEM; repeat=1, ensemblealg=EnsembleThreads(), rng=MersenneTwister[MersenneTwister(1941964948, (0, 15808554, 15807552, 751)), MersenneTwister(3883929895, (0, 15806550, 15805548, 787)), MersenneTwister(5825894842, (0, 13550046, 13549044, 555)), MersenneTwister(7767859789, (0, 13548042, 13547040, 858)), MersenneTwister(9709824736, (0, 13555056, 13554054, 728)), MersenneTwister(11651789683, (0, 13550046, 13549044, 522)), MersenneTwister(13593754630, (0, 13551048, 13550046, 343)), MersenneTwister(15535719577, (0, 13549044, 13548042, 301))])
Closest candidates are:
fit(::Pumas.PumasEMModel, ::AbstractVector{<:Pumas.AbstractSubject}, ::NamedTuple, ::Pumas.SAEM; ensemblealg, diffeq_options, rng) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/bKkN8/src/estimation/saem.jl:3443 got unsupported keyword argument “repeat”
fit(::Pumas.PumasEMModel, ::AbstractVector{<:Pumas.AbstractSubject}, ::NamedTuple, ::Union{Pumas.LikelihoodApproximation, Pumas.MAP}; optimize_fn, constantcoef, ensemblealg, checkidentification, diffeq_options) at /builds/PumasAI/PumasSystemImages-jl/.julia/packages/Pumas/bKkN8/src/estimation/likelihoods.jl:2985 got unsupported keyword arguments “repeat”, “rng”
fit(::Pumas.AbstractPumasModel, ::DataFrame, ::NamedTuple, ::Any…; kwargs…)
All those options, including repeat should be inside Pumas.SAEM()
I see @vijay. It worked. Should I assume the final estimates are geometric mean across repetitions, when repeat>1 ? Or does it mean the estimates of first repeat are used as the initial estimates for the second repeat and so on.
Bob
SAEM is a stochastic algorithm. I think, but I need to confirm, that repeat takes the final estimates of the first repetition and uses them as an initial estimates of the second repetition, and the final estimates of the second are used as an initial of the third, and so on. This way it is a confirmation the the stochastic algorithm has reached a stationary phase. I may be incorrect here, but I will let somebody else comment. @chriselrod
So, to answer your question, these are not the geometric mean estimates of the final parameters
I am getting following error , while I was running SAEM fit :
nested task error: DomainError with -0.0024378399565848086:
sqrt will only return a complex result if called with a complex argument. Try sqrt(Complex(x)).
Please help me to understand the above error.