121
121
# https://github.com/marisae/cancer-chemo-identifiability/blob/master/Profile%20Likelihood/testa0_fit.m#L40-L41
122
122
sigmasq = (mean ([(Cerr005/ C005); (Cerr010/ C010); (Cerr040/ C040); (Cerr100/ C100)]))^ 2
123
123
124
- #=
125
- pid = 2
126
- maxiters = 1000
127
- scan_bounds = ((1 - 0.05*14)*p0[pid], (1 + 0.05*60)*p0[pid])
128
- Δθ = 0.2
129
-
130
- identity_profiler = IntegrationProfiler(integrator=Rodas4(autodiff=false), matrix_type=:identity, integrator_opts=(reltol=1e-3, abstol=1e-6), gamma = 1.)
131
- hess_profiler = IntegrationProfiler(integrator=Rodas4(autodiff=false), matrix_type=:hessian, integrator_opts=(reltol=1e-3, abstol=1e-6), gamma = 1.)
132
-
133
- opt_profiler0 = OptimizationProfiler(NLopt.LN_NELDERMEAD, (reltol=1e-3,), :fixed, Δθ, maxiters)
134
- opt_profiler1 = OptimizationProfiler(NLopt.LN_NELDERMEAD, (reltol=1e-3,), :linear, Δθ, maxiters)
135
-
136
- @time opt_prof0 = profile_parameter(_prob, opt_profiler0, pid; scan_bounds, threashold=threshold, verbose=true)
137
- @time opt_prof1 = profile_parameter(_prob, opt_profiler1, pid; scan_bounds, threashold=threshold, verbose=true)
138
-
139
- @time hess_prof = profile_parameter(_prob, hess_profiler, pid; scan_bounds, threashold=threshold, verbose=true)
140
- @time identity_prof = profile_parameter(_prob, identity_profiler, pid; scan_bounds, threashold=threshold, verbose=true)
141
-
142
- plot(identity_prof, dpi=400, legend=:top, xlabel="ka", ylabel="L(θ)")
143
- plot!(hess_prof, alpha=0, conf_level=false, steps=false)
144
-
145
- plot(opt_prof0, dpi=400, legend = :bottomleft, xlabel="ka", ylabel="L(θ)")
146
- plot(opt_prof1, dpi=400, legend=:bottomleft, xlabel="ka", ylabel="L(θ)")
147
-
148
- endp = LikelihoodProfiler.ProfileSolution(
149
- nothing,
150
- nothing,
151
- hess_prof.θopt,
152
- [[0., 4.92292246374492, 0., 0., 0.], [0., 10.778706398866387, 0., 0., 0.]],
153
- [hess_prof.flevel, hess_prof.flevel],
154
- 2,
155
- #syms::
156
- hess_prof.flevel,
157
- nothing,
158
- nothing
159
-
160
- )
161
-
162
- df = CSV.read("./lp_ka.csv", DataFrame)
163
- rename!(df, :f=>:obj)
164
-
165
- ka_full = LikelihoodProfiler.ProfileSolution(
166
- nothing,
167
- nothing,
168
- hess_prof.θopt,
169
- [[df.ka[i]] for i in 1:64],
170
- [df.obj[i] for i in 1:64],
171
- 1,
172
- #syms::
173
- hess_prof.flevel,
174
- nothing,
175
- nothing
176
-
177
- )
178
-
179
- plot(endp, xticks=[4, 6,7,8,9,10, 4.99, 10.73, 12], xlims=(4.5,11),
180
- dpi=400, legend=:top, xlabel="ka", ylabel="L(θ)", endpoints=true)
181
- scatter!((8.0959,0.538146855), label="optim")
182
- plot!(ka_full, alpha=0.2, conf_level=false)
183
- =#
0 commit comments