diff --git a/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl b/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl index a1bf6fa..7e178cf 100644 --- a/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl +++ b/exercises/solved_notebooks/P5_mcmc/MCMC_2-basics.jl @@ -62,7 +62,7 @@ md"### 2: Unconditional expected value of Y" # ╔═╡ 3decb2ec-210a-4b2f-842d-6fd40dd3f77b @model function mole() - X ~ X_prior + X ~ Exponential(100) Y ~ Uniform(0, X) return Y end @@ -108,7 +108,7 @@ md"### 5: Conditional distribution of X (with more data)" # ╔═╡ b2b16c34-e12a-4bea-8098-313d01913bbf @model function mole2() - X ~ X_prior + X ~ Exponential(100) Ys = zeros(4) for i in 1:length(Ys) Ys[i] ~ Uniform(0, X) @@ -245,7 +245,7 @@ md"### 3: Conditional expected value" # ╔═╡ 126d3954-c77d-4c98-abe0-fd87d14e6265 @model function lights() - μ ~ lights_prior + μ ~ LogNormal(log(40), 0.5) lifespans = zeros(4) for i in 1:length(lifespans) lifespans[i] ~ Exponential(μ) @@ -284,7 +284,7 @@ md""" # ╔═╡ 8fc58fa4-b005-4f32-9eae-a8143582a1ae @model function lights_censored(time_observed, n_working) - μ ~ lights_prior + μ ~ LogNormal(log(40), 0.5) lifespans = zeros(2) for i in 1:length(lifespans) diff --git a/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl b/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl index 60facfa..2c3bdb5 100644 --- a/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl +++ b/exercises/solved_notebooks/P5_mcmc/MCMC_3-advanced.jl @@ -211,12 +211,18 @@ end # ╔═╡ d9b50958-20ae-4085-80d6-19420c7d89df @model function petrigrowth🌟(t_obs) - P0 ~ dropletdist🌟 + P0 ~ MixtureModel( + [ + truncated(Normal(10, sqrt(10)), lower = 0.0), + truncated(Normal(30, sqrt(30)), lower = 0.0) + ], + [0.75, 0.25] + ) r ~ LogNormal(0.0, 0.3) K ~ Normal(1e5, 1e4) logfun = t -> logistic(t, P0, r, K) - Pt = logfun(t_obs) + Pt = max(logfun(t_obs), 0) # set to 0 if negative to prevent possible errors (not required) P_obs ~ Poisson(Pt) return logfun diff --git a/exercises/solved_notebooks/P6_calib/calib_fermenter_monod_sol.jl b/exercises/solved_notebooks/P6_calib/calib_fermenter_monod_sol.jl index 0daed9d..19bb6ba 100644 --- a/exercises/solved_notebooks/P6_calib/calib_fermenter_monod_sol.jl +++ b/exercises/solved_notebooks/P6_calib/calib_fermenter_monod_sol.jl @@ -113,34 +113,23 @@ Furthermore, suppose that at $t = 0\;h$ no substrate $S$ is present in the react Calibrate the parameter values for $\mu_{max}$ and $K_s$ using the aforementioned measurement data for $S$ and $X$ in a timespan of `[0, 100]`$\mathrm{h}$. Take the values above as initial values for $\mu_{max}$ and $K_s$. """ -# ╔═╡ 7a227eaf-18d0-44f4-ac4b-f529e81c7471 -md""" -Create an `ODEProblem`. Use the aforementioned values as initial values for the problem. -""" - -# ╔═╡ 6375478f-1af9-4fd2-b6f3-101a6f796f2d +# ╔═╡ 1689c839-9aa5-4bff-a26f-0c8fbf47aa32 # u0 = missing u0 = [:S=>0.0, :X=>0.0005] -# ╔═╡ 38fe8304-af61-40a7-ac86-480dfb892185 +# ╔═╡ 3fbc29b0-c299-4223-9922-d6fb9a4a6a90 # tspan = missing tspan = (0.0, 100) -# ╔═╡ 87482f88-8413-4820-9613-7941f3d61bd7 -# parms = missing -parms = [:μmax=>0.40, :Ks=>0.015, :Y=>0.67, :Q=>2, :V=>40, :Sin=>0.022] - -# ╔═╡ 94f3bd7b-5c2c-4661-a0ab-2cdaf2cd6743 -# oprob = missing -oprob = ODEProblem(fermenter_monod, u0, tspan, parms, combinatoric_ratelaws=false) - # ╔═╡ f6a8f134-6db0-4d74-8af5-82826347d8f0 md""" -Declare the Turing model. Use `InverseGamma()` for the standard deviations of the measurements and `LogNormal()` for `μmax` and `Ks`. +Declare the Turing model. Assume the following for the priors: +- The measurement error is an unknown positive value, but probably near $0$. +- The parameters `μmax` and `K` are both positive and expected to be in $[0.0, 1.0]$, presumably around $0.1$. """ # ╔═╡ 4c28a66a-ee2c-42a2-95c7-ea4ddb6a232d -# @model function fermenter_inference(t_meas, S, X) +# @model function fermenter_inference(t_meas) # σ_S ~ missing # σ_X ~ missing # μmax ~ missing @@ -151,18 +140,24 @@ Declare the Turing model. Use `InverseGamma()` for the standard deviations of th # S ~ missing # X ~ missing # end -@model function fermenter_inference(t_meas, S, X) - σ_S ~ InverseGamma() - σ_X ~ InverseGamma() - μmax ~ LogNormal() - Ks ~ LogNormal() - parms = [:μmax=>μmax, :Ks=>Ks, :Y=>0.67, :Q=>2, :V=>40, :Sin=>0.022] - oprob = ODEProblem(fermenter_monod, u0, tspan, parms, combinatoric_ratelaws=false) - osol = solve(oprob, Tsit5(), saveat=t_meas) - S ~ MvNormal(osol[:S], σ_S^2 * I) - X ~ MvNormal(osol[:X], σ_X^2 * I) +@model function fermenter_inference(t_meas) + σ_S ~ Exponential() + σ_X ~ Exponential() + μmax ~ LogNormal(log(0.1)) + Ks ~ LogNormal(log(0.1)) + parms = [:μmax => μmax, :Ks => Ks, :Y => 0.67, :Q => 2, :V => 40, :Sin => 0.022] + oprob = ODEProblem(fermenter_monod, u0, tspan, parms) + osol = solve(oprob, AutoTsit5(Rosenbrock23()), saveat=t_meas) + S ~ MvNormal(osol[:S], σ_S) + X ~ MvNormal(osol[:X], σ_X) end +# ╔═╡ c2120cce-5cae-42e5-b1c6-1d10b49d9ffc +md""" +!!! tip + This model can change between being non-stiff and stiff based on the sampled parameter values. You can use an auto-switching solver such as `AutoTsit5(Rosenbrock23())` here to make calibration more stable. +""" + # ╔═╡ 3136b15d-5078-4bcd-954b-e89bcb8aed1b md""" Provide the measurements to the Turing model. @@ -170,15 +165,15 @@ Provide the measurements to the Turing model. # ╔═╡ 6a508a62-61b9-4273-8e45-b26f594e8da9 # fermenter_inf = missing -fermenter_inf = fermenter_inference(t_meas, S_meas, X_meas) +fermenter_inf = fermenter_inference(t_meas) | (S = S_meas, X = X_meas) # ╔═╡ 63420055-55f8-4def-8b0e-11ea61483010 md""" -Optimize the priors ($\sigma_S$, $\sigma_X$, $\mu_{max}$ and $K_s$). Do this with `MLE` method and Nelder-Mead. Store the optimization results in `results_mle`. +Optimize the likelihood of the parameters ($\sigma_S$, $\sigma_X$, $\mu_{max}$ and $K_s$) using the NelderMead optimizer. Store the optimization results in `results_mle`. """ # ╔═╡ d52c9da8-d8a4-4db0-ac6d-6d16ccf4775c -# results_map = missing +# results_mle = missing # Uncomment and complete the instruction results_mle = optimize(fermenter_inf, MLE(), NelderMead()) # ╔═╡ e1b0ee01-f16c-40e9-a0f9-80072d690936 @@ -223,8 +218,8 @@ Create an ODEProblem and solve it. Use `Tsit5()` and `saveat=0.5`. """ # ╔═╡ 853c1a92-d50f-4b05-9ed3-d3ee1656665a -# oprob_opt = missing -oprob_opt = ODEProblem(fermenter_monod, u0, tspan, parms_opt, combinatoric_ratelaws=false) +# oprob_opt = missing # Uncomment and complete the instruction +oprob_opt = ODEProblem(fermenter_monod, u0, tspan, parms_opt) # ╔═╡ f45e8124-e942-438e-99c5-3032ccc01454 # osol_opt = missing @@ -266,13 +261,11 @@ end # ╟─6c481447-28c6-4530-bf2c-64762121bc71 # ╠═918fd524-81fa-4aff-a403-37402e47235b # ╟─ef977370-06ee-4a73-85e2-609a744167d3 -# ╟─7a227eaf-18d0-44f4-ac4b-f529e81c7471 -# ╠═6375478f-1af9-4fd2-b6f3-101a6f796f2d -# ╠═38fe8304-af61-40a7-ac86-480dfb892185 -# ╠═87482f88-8413-4820-9613-7941f3d61bd7 -# ╠═94f3bd7b-5c2c-4661-a0ab-2cdaf2cd6743 +# ╠═1689c839-9aa5-4bff-a26f-0c8fbf47aa32 +# ╠═3fbc29b0-c299-4223-9922-d6fb9a4a6a90 # ╟─f6a8f134-6db0-4d74-8af5-82826347d8f0 # ╠═4c28a66a-ee2c-42a2-95c7-ea4ddb6a232d +# ╟─c2120cce-5cae-42e5-b1c6-1d10b49d9ffc # ╟─3136b15d-5078-4bcd-954b-e89bcb8aed1b # ╠═6a508a62-61b9-4273-8e45-b26f594e8da9 # ╟─63420055-55f8-4def-8b0e-11ea61483010 diff --git a/exercises/solved_notebooks/P6_calib/calib_intro.jl b/exercises/solved_notebooks/P6_calib/calib_intro.jl index 1cfea7a..bdefe5f 100644 --- a/exercises/solved_notebooks/P6_calib/calib_intro.jl +++ b/exercises/solved_notebooks/P6_calib/calib_intro.jl @@ -4,6 +4,18 @@ using Markdown using InteractiveUtils +# This Pluto notebook uses @bind for interactivity. When running this notebook outside of Pluto, the following 'mock version' of @bind gives bound variables a default value (instead of an error). +macro bind(def, element) + #! format: off + return quote + local iv = try Base.loaded_modules[Base.PkgId(Base.UUID("6e696c72-6542-2067-7265-42206c756150"), "AbstractPlutoDingetjes")].Bonds.initial_value catch; b -> missing; end + local el = $(esc(element)) + global $(esc(def)) = Core.applicable(Base.get, el) ? Base.get(el) : iv(el) + el + end + #! format: on +end + # ╔═╡ f8a92690-990b-4341-89e1-322adbcb8d1b using Pkg; Pkg.activate("..") @@ -42,18 +54,22 @@ In the models discussed in the previous sessions, we always knew the values of a # ╔═╡ 987f0a4d-e416-4ceb-adbe-3dcdca9d0996 md""" -The search of optimal parameter values usually involves a function, such as a loss function, a (log) likelihood function or a posterior distribution function. In this session we will be (mainly) using the MLE (Maximum Likelihood Estimation) and MAP (Maximum A Posteriori estimation) methods. +The search of optimal parameter values usually involves a function, such as a loss function, a (log) likelihood function or a posterior distribution function. In this session we will be (mainly) using the MLE (Maximum Likelihood Estimation) and MAP (Maximum A Posteriori estimation) methods, which respectively try to maximize the likelihood function and the posterior probability of the data. -In the MLE method, a likelihood function of a given probability density function is maximized during the search of optimal parameter values of a model in order to fit experimental data. The parameter values are considered unknown but viewed as fixed points. +In order to understand the difference, recall Bayes' theorem applied to a set of parameters $θ$ and data $D$: +```math +P(θ \mid D) = \frac{P(D \mid θ) \, P(θ)}{P(D)} +``` +In the MLE method, the likelihood function $P(D \mid θ)$ (= the probability of the data given the parameters) is maximized during the search of optimal parameter values of a model in order to fit experimental data. The parameter values are considered unknown but viewed as fixed points. -In the MAP method, a posterior distribution function is maximized. Instead of viewing the parameter values as fixed points, they are now treated as random variables in the model which follow a prior distribution. In other words, we have prior belief in which distribution these parameters come from (Normal, Beta, etc). Once new data comes in, we update our prior belief, leading to a posterior belief. Hence, we now have a better idea from which distribution these parameters come. +In the MAP method, the posterior probability $P(θ \mid D)$ (= the probability of the parameters given the data) is maximized instead. Instead of viewing the parameter values as fixed points, they are now treated as random variables in the model which follow a prior distribution. In other words, we have prior belief in which distribution these parameters come from (Normal, Beta, etc). Once new data comes in, we update our prior belief, leading to a posterior belief. Hence, we now have a better idea from which distribution these parameters come. -One caveat with the MAP method is that it only considers the most likely point without taking into account other values of parameters from posterior distribution, which leads to a huge loss of useful information in posterior distribution. A better, yet computationally exhaustive method is using the MCMC (Markov chain Monte Carlo) sampling methods. Here we will use the NUTS sampler in this session. +One caveat with the MLE and MAP methods are that they only return a point estimate of the optimal value, which gives no information on how certain we are about this value. A more informative, yet computationally exhaustive method is using the MCMC (Markov chain Monte Carlo) sampling methods, which approximate the entire posterior distribution of the unknown values. """ # ╔═╡ 75efff36-8da7-4d04-afa2-a2f8324bc103 md""" -In this notebook we will calibrate the different parameters involved in the grass growth models. Therefore, we will implement an objective function for each model and then minimizing them using an optimization algorithm. To illustrate this concept, we first revisit the three simple models modelling the grass growth yield. +In this notebook we will calibrate the different parameters involved in the grass growth models. To illustrate this concept, we first revisit the three simple models modelling the grass growth yield. """ # ╔═╡ 3dcb9c9d-370b-4031-b7c0-cee80742557a @@ -69,7 +85,7 @@ In this notebook, three different models will be used, each modelling the yield - Exponential growth model: $\cfrac{dW}{dt} = \mu \left( W_f - W \right)$ - Gompertz growth model: $\cfrac{dW}{dt} = \left( \mu - d \ln(W) \right) W$ -with output $W$ the grass yield, and $W_f$, $\mu$ and $d$ parameters. The table below show the parameter values and the initial condition that will be used as initial values for the optimization algorithm: +with output $W$ the grass yield, and $W_f$, $\mu$ and $d$ parameters. The table below shows some typical parameter values and initial conditions for grasslands similar to the one we observed, which we can use as prior information. | | $\mu$ | $W_f$ | $d$ | $W_0$ | |:----------- |:----------:|:-----------:|:------------:|:------------:| @@ -93,133 +109,180 @@ md""" Variables containing the initial condition and parameters values will be defined later in the objective function. """ -# ╔═╡ 9a5bc72b-346d-4e95-a873-783037ed98bc +# ╔═╡ 1aa44f2b-6f33-437f-b9dd-89762d9f28ea md""" -### Logistic growth model +### The measurement data """ -# ╔═╡ 4b17f7e6-26b0-425a-9f53-ccf1689639fe +# ╔═╡ b2b433ed-0266-4bea-a7e8-32adba542d4c md""" -We will illustrate the calibration with the logistic growth model. The latter can be done via a Catalyst reaction network model or via a model built with ModelingToolkit. As an example, we will show both possibilities for the logistic growth model. +Assume that the measured grass yields (of a certain plant type) are the following: """ -# ╔═╡ 6736542c-5378-480e-a56b-65956b416225 +# ╔═╡ 7c966a66-0091-4b81-9a7e-02ccd0d3db10 +W_meas = [1.87, 2.45, 3.72, 4.32, 5.28, 7.01, 6.83, 8.62, 9.45, 10.31, 10.56, 11.72, 11.05, 11.53, 11.39, 11.7, 11.15, 11.49, 12.04, 11.95, 11.68] + +# ╔═╡ 3edd2acc-a865-4675-afef-8868c68256f1 md""" -#### 1) Catalyst based model +They have been measured at the following corresponding time instances: """ -# ╔═╡ ba56adb1-9405-40d5-be48-4273b42ab145 -# growth_log = @reaction_network begin -# @species W(t)=2.0 -# @parameters μ=0.07 Wf=10.0 -# μ*(1-W/Wf), W --> 2W -# end +# ╔═╡ 877298e8-b61b-4c3a-ba2c-2827acdcfb50 +t_meas = 0:5:100 -# ╔═╡ 1c311b6f-6b18-4170-b500-33a8e4d3cb2d +# ╔═╡ ef06cc43-510b-4ff9-b0b7-1c7fc267e9b1 md""" -#### ̇2) ModelingToolkit based model +We can make a scatter plot of this data (including a title, a legend label, an X-axis label, X- and Y-axis limits) in the following way: """ -# ╔═╡ 50220e1f-8e42-44da-be03-08c11df967f0 -@variables W(t) - -# ╔═╡ 449286b0-0210-4c6b-b28f-fca87b52d674 -@parameters μ Wf d +# ╔═╡ cb2bc6ee-4211-47e1-9956-5cf1b0c0671d +scatter(t_meas, W_meas, title="Grass growth data", + label="Yield", + xlabel="t", + xlims=(0, 100), + ylims=(0, 14)) -# ╔═╡ f720eac9-cb29-4eef-ac94-471395941c0f -# IF YOU UNCOMMENT THE FOLLOWING LINE, PLEASE COMMENT THE CATALYST MODEL FIRST!!! -@mtkbuild growth_log = ODESystem([D(W) ~ μ*(1 - W/Wf)*W], t) +# ╔═╡ 50e6af5b-04f4-495c-853b-c746fd27254f +md"## Calibration with Turing" -# ╔═╡ 024cce41-0469-415f-bdc5-69e8d38c4b69 +# ╔═╡ a665740a-ea94-4452-8cbc-a1647319dfce md""" -#### Creating the ODEProblem +We will be using the familiar Turing framework to perform optimisation, which consists roughly of the following steps: +- Define your model as a Turing model, with the following components: + - The observed input data as input of the model + - Priors for all unknown values such as initial values, parameter values and the measurement error + - An (ODE-based) model that predicts values of the output variable based on the observed input values and the unknown initial values, parameter values, etc. + - The relationship between the observed outputs and the predicted outputs, which often comes down to specifying the measurement error +- Instantiate the Turing model and condition it on the observed values of the output variable +- Call an optimisation algorithm on the Turing model +- Extract the optimized parameters +- Visualise the results """ -# ╔═╡ acccb2fa-12b2-4fc7-91e3-58a4b1a02892 +# ╔═╡ b2f20a9a-0613-4bb4-a1d3-df1620aef105 +md"### The priors" + +# ╔═╡ 6852f166-350e-4434-8304-382c01113ba1 md""" -Next, we will need to create an `ODEProblem` in advance before we can optimize some of its parameters. We will provide the values in the aforementioned table as initial values for the problem. +During calibration, the prior distributions have two roles: +1. Define for every unknown quantity **the bounds**: the optimizer will only search within the domain of the prior distribution. Note: you can add additional bounds to any existing distribution using the `truncated` function. +1. Define part of the **posterior probability** of observing the data given the parameter values. In other words, the theoretically optimal values depend on your prior distributions. + +**Note that the second only applies to methods that take prior information into account, which Maximum Likelihood Estimation (MLE) does not.** """ -# ╔═╡ e54ea8d1-0854-44fa-aed8-45d106e921e4 -u0_log = [:W=>2.0] +# ╔═╡ 4f78b2e8-9132-4175-b9ae-9f5979649e90 +md"#### Choice of priors for the grass growth models" -# ╔═╡ 8d96eb17-ce19-4523-916f-3cd0441a16ca -parms_log = [:μ=>0.07, :Wf=>10.0] +# ╔═╡ 557e69c5-d0c6-41ee-b3df-851bd1d16923 +md""" +In general for our grass models, we will need to define priors for the following: +- the measurement error (standard deviation) $\sigma_W$. +- the initial condition $W_0$. +- the parameters $\mu$ and either $W_f$ or $d$ (depending on the model). +""" -# ╔═╡ 561f97d4-b7c3-491a-9031-5148336867ea +# ╔═╡ 237934b7-c89a-4ef6-938f-4629c0fa4c95 md""" -Now you can choose which model (based on Catalyst or MTK) to use when creating the ODEProblem. They should both render the same results. *Leave one of them commented!* +For the measurement error $\sigma_W$, a distribution with most of its probability density around 0 and a long positive tail is often a good choice. We will use the Exponential distribution here for this purpose. As for its parameters, considering our yield is in the order of 1 to 10 [t/ha], an expected value (and therefore standard deviation) of around 1 seems appropriate. """ -# ╔═╡ 5c9db9df-0cbd-41ac-afe9-fb5616c967be -oprob_log = ODEProblem(growth_log, u0_log, tspan, parms_log); +# ╔═╡ 513ee94d-5d49-43d2-a1db-638fe9e2de1a +@bind example_mean_error Slider(0.01:0.01:10.0, show_value = true, default = 1) -# ╔═╡ 1aa44f2b-6f33-437f-b9dd-89762d9f28ea +# ╔═╡ 8e163f77-b762-427a-b005-313284dc2e9e +plot(Exponential(example_mean_error), title = "Example prior for the measurement error", legend = false) + +# ╔═╡ d9c38e9d-12b8-4e98-8546-bbc123960e6f md""" -### The measurement data +For the initial- and parameter values, we will default to log-normal distributions. These have the following properties that generally fit well with biological parameters, such as is the case in this notebook: +- A positive domain, which bounds possible values to the positive numbers +- Most of the probability density around the expected value +- A long postive tail, which allows for outliers + +As the distribution of a variable whose logarithm is normally distributed, it does have some strange behaviour: the arguments of the distribution are the mean and standard deviation of the **logarithm**. To clarify, for $X \sim \mathrm{LogNormal(μ, σ)}$: +- $μ = E[\mathrm{log}(X)]$ +- $σ = \sqrt{\mathrm{Var}(\mathrm{log}(X))}$ + +Practically, this means **you need to specify the logarithm of the desired expected value and play around with the standard deviation**. For this exercise, we will choose an expected value based on the table discussed earlier in this section and keep the standard deviation to the default value of 1. Additionally, for some parameters we may want to truncate them to to keep them within (biologically and numerically) sensible bounds. """ -# ╔═╡ b2b433ed-0266-4bea-a7e8-32adba542d4c +# ╔═╡ 39c70e22-86e0-467e-abab-24f99548b096 +@bind example_log_μ Slider(0.01:0.01:10.0, default = 0.1, show_value = true) + +# ╔═╡ ddef987d-59db-4adf-ae77-d4d32a41ddec +@bind example_log_σ Slider(0.01:0.01:3.0, default = 1.0, show_value = true) + +# ╔═╡ 54fb706b-9036-42b5-833e-4081f4c2179c +begin + example_param_prior = LogNormal(log(example_log_μ), example_log_σ) + + println("""Look at this strange distribution with + - mean = $(mean(example_param_prior)) + - median = $(median(example_param_prior)) + - σ = $(std(example_param_prior))""") + plot(example_param_prior, xlims = (0, 5*example_log_μ), title = "Example prior for the model parameters", legend = false) +end + +# ╔═╡ 9a5bc72b-346d-4e95-a873-783037ed98bc md""" -Assume that the measured grass yields (of a certain plant type) are the following: +## Example: the Logistic growth model """ -# ╔═╡ 7c966a66-0091-4b81-9a7e-02ccd0d3db10 -W_meas = [1.87, 2.45, 3.72, 4.32, 5.28, 7.01, 6.83, 8.62, 9.45, 10.31, 10.56, 11.72, 11.05, 11.53, 11.39, 11.7, 11.15, 11.49, 12.04, 11.95, 11.68] +# ╔═╡ 4b17f7e6-26b0-425a-9f53-ccf1689639fe +md""" +We will illustrate the calibration with the logistic growth model. The latter can be done via a Catalyst reaction network model or via a model built with ModelingToolkit. As an example, we will show both possibilities for the logistic growth model. +""" -# ╔═╡ 3edd2acc-a865-4675-afef-8868c68256f1 +# ╔═╡ 6736542c-5378-480e-a56b-65956b416225 md""" -They have been measured at the following corresponding time instances: +#### 1) Catalyst based model """ -# ╔═╡ 877298e8-b61b-4c3a-ba2c-2827acdcfb50 -t_meas = 0:5:100 +# ╔═╡ ba56adb1-9405-40d5-be48-4273b42ab145 +# growth_log = @reaction_network begin +# @species W(t)=2.0 +# @parameters μ=0.07 Wf=10.0 +# μ*(1-W/Wf), W --> 2W +# end -# ╔═╡ ef06cc43-510b-4ff9-b0b7-1c7fc267e9b1 +# ╔═╡ 1c311b6f-6b18-4170-b500-33a8e4d3cb2d md""" -We can make a scatter plot of this data (including a title, a legend label, an X-axis label, X- and Y-axis limits) in the following way: +#### ̇2) ModelingToolkit based model """ -# ╔═╡ cb2bc6ee-4211-47e1-9956-5cf1b0c0671d -scatter(t_meas, W_meas, title="Grass growth data", - label="Yield", - xlabel="t", - xlims=(0, 100), - ylims=(0, 14)) +# ╔═╡ 50220e1f-8e42-44da-be03-08c11df967f0 +@variables W(t) + +# ╔═╡ 449286b0-0210-4c6b-b28f-fca87b52d674 +@parameters μ Wf d + +# ╔═╡ f720eac9-cb29-4eef-ac94-471395941c0f +# IF YOU UNCOMMENT THE FOLLOWING LINE, PLEASE COMMENT THE CATALYST MODEL FIRST!!! +@mtkbuild growth_log = ODESystem([D(W) ~ μ*(1 - W/Wf)*W], t) # ╔═╡ d75246d4-e03b-4684-be7d-4bcfb61ed7ef md""" ### Declaration of the Turing model """ -# ╔═╡ dd3a32f1-bdb6-44a3-acbe-f4269725c9e4 -md""" -In the Turing model we will define our priors for the following magnitudes: -- the measurement error (standard deviation) $\sigma_W$, -- the initial condition $W_0$, and -- the parameters $\mu$ and $W_f$. - -We will thereby take an Inverse Gamma prior distribution for $\sigma_W$ and LogNormal prior distributions for the initial condition and the parameters (assuming they lay in pre-defined range). -""" - # ╔═╡ 8a9115eb-4044-4cab-a7db-39b5dd86c70d -@model function growth_log_inference(t_meas, W_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - Wf ~ LogNormal() +@model function growth_log_inference(t_meas) + σ_W ~ Exponential() + W0 ~ LogNormal(log(1)) + μ ~ LogNormal(log(0.1)) + Wf ~ truncated(LogNormal(log(10)), lower = 0.1) # prevent zeros in denominator of our model u0_log = [:W => W0] parms_log = [:μ => μ, :Wf => Wf] oprob_log = ODEProblem(growth_log, u0_log, tspan, parms_log) osol_log = solve(oprob_log, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_meas ~ MvNormal(osol_log[:W], σ_W^2 * I) + W_s ~ MvNormal(osol_log[:W], σ_W) end # ╔═╡ 48c9f616-d298-40da-b917-225abd39b3d9 md""" Some remarks: - The time points are the ones from the measurements, therefore, we set: `saveat=t_meas`. -- We need to solve the ODE problem inside the Turing model function with values for $W_0$, $\mu$ and $W_f$ sampled from the distributions. Therefore, we need to remake our ODE problem with the appropriate initial and parameter values and solve it. - Depending on the priors, the parameter values of our ODE may vary wildly during calibration. As this can influence the stiffness of the system, it can be beneficial to use an ODE solver that automatically detects the stiffness of the system and switches to a stiff solver if necessary, such as `AutoTsit5(Rosenbrock23())`. We don't expect you to know when this is necessary, but if your calibration regularly gives instability errors, this may help! """ @@ -229,7 +292,7 @@ We will provide the measurements to the Turing model: """ # ╔═╡ 8b6534d6-776b-4285-8498-a9b34051facc -growth_log_inf = growth_log_inference(t_meas, W_meas) +growth_log_inf = growth_log_inference(t_meas) | (W_s = W_meas, ) # ╔═╡ a6972aef-63ad-401c-acf5-6d59f9fc6698 md""" @@ -407,17 +470,20 @@ We will use Markov chain Monte Carlo (MCMC) method in combination with the No U- # ╔═╡ 0c047043-3284-422a-9c88-2f4f4c170edf results_log_nuts = sample(growth_log_inf, NUTS(), 500) +# ╔═╡ 78a608a4-0ba7-400e-abb9-a9285c60681c +plot(results_log_nuts) # check convergence + # ╔═╡ 19c362cb-2764-41c9-a571-2e8e2bfcde93 # summarize(results_log_nuts) # ╔═╡ 93db47b2-34e8-43b4-beac-b5620fd444e7 -W0_opt3_log = summarize(results_log_nuts)[:W0, :mean] +W0_opt3_log = mean(results_log_nuts[:W0]) # ╔═╡ 68c71cd2-6cdc-4b80-b6e7-75bfea344295 -μ_opt3_log = summarize(results_log_nuts)[:μ ,:mean] +μ_opt3_log = mean(results_log_nuts[:μ]) # ╔═╡ d6b9eaba-1d43-4e56-8fa1-bb2f87f6fe79 -Wf_opt3_log = summarize(results_log_nuts)[:Wf ,:mean] +Wf_opt3_log = mean(results_log_nuts[:Wf]) # ╔═╡ b843ee5e-9618-4b50-93db-77e19b4be366 md""" @@ -499,31 +565,9 @@ Convert the reaction network model into an ODE system to verify. # missing convert(ODESystem, growth_exp) -# ╔═╡ c6d373f4-f13c-4135-823d-ee8fbeb71b56 -md""" -Create an `ODEProblem`. Use the values in the aforementioned table as initial values for the problem. Use the same `tspan` as before. -""" - -# ╔═╡ a97abaa7-b642-4201-86f1-5c8995b07536 -# u₀_exp = missing -u0_exp = [:W=>2.0] - -# ╔═╡ 387730b4-bd06-492f-94e6-231bd68b3436 -# parms_exp = missing -parms_exp = [:μ=>0.02, :Wf=>10.0] - -# ╔═╡ 290a7fe8-3b1e-423f-8b30-9bd8903d2e8f -# oprob_exp = missing -oprob_exp = ODEProblem(growth_exp, u0_exp, tspan, parms_exp) - -# ╔═╡ febe2b67-2a8f-4575-946d-30877bd5f2d4 -md""" -Use the same measurement data (`W_meas`, `t_meas`) as before. -""" - # ╔═╡ b4300e8a-8052-419b-98c8-0508ebee2393 md""" -Declare the Turing model. Take the same priors (and distributions) as before. +Declare the Turing model. Try to find sensible prior distributions. """ # ╔═╡ 2c6ae74c-2da4-4867-ad8a-f4e835101d63 @@ -538,16 +582,16 @@ Declare the Turing model. Take the same priors (and distributions) as before. # osol_exp = missing # W_meas ~ missing # end -@model function growth_exp_inference(t_meas, W_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - Wf ~ LogNormal() +@model function growth_exp_inference(t_meas) + σ_W ~ Exponential() + W0 ~ LogNormal(log(1)) + μ ~ LogNormal(log(0.1)) + Wf ~ LogNormal(log(10)) u0_exp = [:W => W0] parms_exp = [:μ => μ, :Wf => Wf] oprob_exp = ODEProblem(growth_exp, u0_exp, tspan, parms_exp) osol_exp = solve(oprob_exp, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_meas ~ MvNormal(osol_exp[:W], σ_W^2 * I) + W_s ~ MvNormal(osol_exp[:W], σ_W) end # ╔═╡ c6a5d453-d610-4f65-847c-c878dd41726c @@ -557,7 +601,7 @@ Provide the measurements to the Turing model. # ╔═╡ fff17cf7-173d-4f64-94a9-4bf46acc882d # growth_exp_inf = missing -growth_exp_inf = growth_exp_inference(t_meas, W_meas) +growth_exp_inf = growth_exp_inference(t_meas) | (W_s = W_meas, ) # ╔═╡ eee55784-a641-445e-be75-0b19e2a94754 md""" @@ -662,33 +706,10 @@ Implement a system using MTK for the Gompertz growth model.\ # ╔═╡ 6e85c08e-9a97-4897-9c81-89517f55e254 @mtkbuild growth_gom = ODESystem([D(W) ~ (μ - d*log(W))*W], t) -# ╔═╡ 47fb9e4c-df6a-4811-9980-99d595a34908 -md""" -Create an `ODEProblem`. Use the values in the aforementioned table as initial values for the problem. Use the same `tspan` as before. -""" - -# ╔═╡ bbd150de-ff9a-4127-a0dd-2f9762f92b07 -# u0_gom = missing -u0_gom = [:W=>2.0] - -# ╔═╡ da5a0cbb-b033-46f1-a300-3954de138835 -# parms_gom = missing -parms_gom = [:μ=>0.09, :d=>0.04] - -# ╔═╡ 73da8f53-c3af-43b0-9b23-60471f1e3587 -# oprob_gom = missing -oprob_gom = ODEProblem(growth_gom, u0_gom, tspan, parms_gom) - -# ╔═╡ b0e67564-efe8-4fb2-bcf2-a711b770244e -md""" -Use the same measurement data (`W_meas`, `t_meas`) as before. -""" - # ╔═╡ e5081280-d226-4834-8932-c89becd8313c md""" -Declare the Turing model. Take the same priors as before. +Declare the Turing model. Try to find sensible prior distributions. """ -# Take for $\sigma_W$ and $W_0$ the same priors (and distributions) as before, but take for $\mu$ a Uniform prior distribution in the range $[0, 2]$ and the same for $D$ but in the range $[0, 1]$. # ╔═╡ c739a908-2353-4e7a-8fbd-f640dc8cabe0 # @model function growth_gom_inference(t_meas, W_meas) @@ -702,16 +723,16 @@ Declare the Turing model. Take the same priors as before. # osol_gom = missing # W_meas ~ missing # end -@model function growth_gom_inference(t_meas, W_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - d ~ LogNormal() +@model function growth_gom_inference(t_meas) + σ_W ~ Exponential() + W0 ~ truncated(LogNormal(log(1.0)), lower = 1e-5) # prevent log(0) + μ ~ truncated(LogNormal(log(0.1)), upper = 1.0) # prevent overly large exponential growth from crashing the solver - this model has no carrying capacity to slow things down + d ~ LogNormal(log(0.1)) u0_gom = [:W => W0] parms_gom = [:μ => μ, :d => d] oprob_gom = ODEProblem(growth_gom, u0_gom, tspan, parms_gom) osol_gom = solve(oprob_gom, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_meas ~ MvNormal(osol_gom[:W], σ_W^2 * I) + W_s ~ MvNormal(osol_gom[:W], σ_W) end # ╔═╡ 1d0383ad-54d6-4ff2-8555-def83bfff0e6 @@ -721,7 +742,7 @@ Provide the measurements to the Turing model. # ╔═╡ cd1cf2f8-9f7f-4ed4-9cb7-1a6efee68ab4 # growth_gom_inf = missing -growth_gom_inf = growth_gom_inference(t_meas, W_meas) +growth_gom_inf = growth_gom_inference(t_meas) | (W_s = W_meas, ) # ╔═╡ aba74ee0-0163-4e15-8b49-d8dcad4839f7 md""" @@ -836,6 +857,26 @@ md"- Answer: missing" # ╟─85cd60a8-b448-4375-9b6d-399c4336c319 # ╠═5b320989-3e0b-447b-bc9a-25fb221ce609 # ╟─2481cd4f-0efc-4450-ab3d-4a5492597f36 +# ╟─1aa44f2b-6f33-437f-b9dd-89762d9f28ea +# ╟─b2b433ed-0266-4bea-a7e8-32adba542d4c +# ╠═7c966a66-0091-4b81-9a7e-02ccd0d3db10 +# ╟─3edd2acc-a865-4675-afef-8868c68256f1 +# ╠═877298e8-b61b-4c3a-ba2c-2827acdcfb50 +# ╟─ef06cc43-510b-4ff9-b0b7-1c7fc267e9b1 +# ╠═cb2bc6ee-4211-47e1-9956-5cf1b0c0671d +# ╟─50e6af5b-04f4-495c-853b-c746fd27254f +# ╟─a665740a-ea94-4452-8cbc-a1647319dfce +# ╟─b2f20a9a-0613-4bb4-a1d3-df1620aef105 +# ╟─6852f166-350e-4434-8304-382c01113ba1 +# ╟─4f78b2e8-9132-4175-b9ae-9f5979649e90 +# ╟─557e69c5-d0c6-41ee-b3df-851bd1d16923 +# ╟─237934b7-c89a-4ef6-938f-4629c0fa4c95 +# ╠═513ee94d-5d49-43d2-a1db-638fe9e2de1a +# ╠═8e163f77-b762-427a-b005-313284dc2e9e +# ╟─d9c38e9d-12b8-4e98-8546-bbc123960e6f +# ╠═39c70e22-86e0-467e-abab-24f99548b096 +# ╠═ddef987d-59db-4adf-ae77-d4d32a41ddec +# ╠═54fb706b-9036-42b5-833e-4081f4c2179c # ╟─9a5bc72b-346d-4e95-a873-783037ed98bc # ╟─4b17f7e6-26b0-425a-9f53-ccf1689639fe # ╟─6736542c-5378-480e-a56b-65956b416225 @@ -844,21 +885,7 @@ md"- Answer: missing" # ╠═50220e1f-8e42-44da-be03-08c11df967f0 # ╠═449286b0-0210-4c6b-b28f-fca87b52d674 # ╠═f720eac9-cb29-4eef-ac94-471395941c0f -# ╟─024cce41-0469-415f-bdc5-69e8d38c4b69 -# ╟─acccb2fa-12b2-4fc7-91e3-58a4b1a02892 -# ╠═e54ea8d1-0854-44fa-aed8-45d106e921e4 -# ╠═8d96eb17-ce19-4523-916f-3cd0441a16ca -# ╟─561f97d4-b7c3-491a-9031-5148336867ea -# ╠═5c9db9df-0cbd-41ac-afe9-fb5616c967be -# ╟─1aa44f2b-6f33-437f-b9dd-89762d9f28ea -# ╟─b2b433ed-0266-4bea-a7e8-32adba542d4c -# ╠═7c966a66-0091-4b81-9a7e-02ccd0d3db10 -# ╟─3edd2acc-a865-4675-afef-8868c68256f1 -# ╠═877298e8-b61b-4c3a-ba2c-2827acdcfb50 -# ╟─ef06cc43-510b-4ff9-b0b7-1c7fc267e9b1 -# ╠═cb2bc6ee-4211-47e1-9956-5cf1b0c0671d # ╟─d75246d4-e03b-4684-be7d-4bcfb61ed7ef -# ╟─dd3a32f1-bdb6-44a3-acbe-f4269725c9e4 # ╠═8a9115eb-4044-4cab-a7db-39b5dd86c70d # ╟─48c9f616-d298-40da-b917-225abd39b3d9 # ╟─35f158c1-858d-4e4d-ac3d-bf4807dad9a0 @@ -905,6 +932,7 @@ md"- Answer: missing" # ╟─29170e2a-9916-438e-92ca-9f4783397b5e # ╟─7ff9fe52-156b-4a92-9058-781670de3abb # ╠═0c047043-3284-422a-9c88-2f4f4c170edf +# ╠═78a608a4-0ba7-400e-abb9-a9285c60681c # ╠═19c362cb-2764-41c9-a571-2e8e2bfcde93 # ╠═93db47b2-34e8-43b4-beac-b5620fd444e7 # ╠═68c71cd2-6cdc-4b80-b6e7-75bfea344295 @@ -925,11 +953,6 @@ md"- Answer: missing" # ╠═cf1a144e-09e9-42a3-b2a3-b8676a200a39 # ╟─d69ed5dd-80cd-427e-9245-e42576a4688d # ╠═8b429656-1c57-4278-9633-71311f400ed8 -# ╟─c6d373f4-f13c-4135-823d-ee8fbeb71b56 -# ╠═a97abaa7-b642-4201-86f1-5c8995b07536 -# ╠═387730b4-bd06-492f-94e6-231bd68b3436 -# ╠═290a7fe8-3b1e-423f-8b30-9bd8903d2e8f -# ╟─febe2b67-2a8f-4575-946d-30877bd5f2d4 # ╟─b4300e8a-8052-419b-98c8-0508ebee2393 # ╠═2c6ae74c-2da4-4867-ad8a-f4e835101d63 # ╟─c6a5d453-d610-4f65-847c-c878dd41726c @@ -955,11 +978,6 @@ md"- Answer: missing" # ╟─785d500b-f8ea-446a-9952-2a5fd5d83d24 # ╟─e754826a-7411-4072-b0dc-a4bad7a15f98 # ╠═6e85c08e-9a97-4897-9c81-89517f55e254 -# ╟─47fb9e4c-df6a-4811-9980-99d595a34908 -# ╠═bbd150de-ff9a-4127-a0dd-2f9762f92b07 -# ╠═da5a0cbb-b033-46f1-a300-3954de138835 -# ╠═73da8f53-c3af-43b0-9b23-60471f1e3587 -# ╟─b0e67564-efe8-4fb2-bcf2-a711b770244e # ╟─e5081280-d226-4834-8932-c89becd8313c # ╠═c739a908-2353-4e7a-8fbd-f640dc8cabe0 # ╟─1d0383ad-54d6-4ff2-8555-def83bfff0e6 diff --git a/exercises/solved_notebooks/P6_calib/calib_irrigation_sol.jl b/exercises/solved_notebooks/P6_calib/calib_irrigation_sol.jl index 18cf1a2..0eab9d9 100644 --- a/exercises/solved_notebooks/P6_calib/calib_irrigation_sol.jl +++ b/exercises/solved_notebooks/P6_calib/calib_irrigation_sol.jl @@ -110,13 +110,13 @@ We can make a scatter plot of the measured data for both $S_1$ and $S_2$ for the # ╔═╡ cef4b9a8-b5bf-4a2d-8a8f-5d8f85534859 begin - scatter(t_meas, S1_meas1, label="S1 meas", color=:blue, title="Experment 1") + scatter(t_meas, S1_meas1, label="S1 meas", color=:blue, title="Experiment 1") scatter!(t_meas, S2_meas1, label="S2 meas", color=:red, ylims=(0, 150)) end # ╔═╡ fc2cabd7-e778-4211-bf87-b5c11ca054c9 begin - scatter(t_meas, S1_meas2, label="S1 meas", color=:blue, title="Experment 2") + scatter(t_meas, S1_meas2, label="S1 meas", color=:blue, title="Experiment 2") scatter!(t_meas, S2_meas2, label="S2 meas", color=:red, ylims=(0, 150)) end @@ -125,34 +125,17 @@ md""" Calibrate the parameter values for $k$ and $S_{max}$ using the aforementioned measurement data for $S_1$ and $S_2$ in a timespan of $[0, 150]\,h$. Take the values from above as initial values for $k$ and $S_{max}$. """ -# ╔═╡ d35bbe54-5ebc-4ea5-a6c8-a6419476ec4c -md""" -Create an `ODEProblem`. Use the aforementioned values as initial values for the problem. -""" - -# ╔═╡ 10057510-e4b6-4a3e-9d3f-f05effc88a58 -# u0 = missing -u0 = [:S₁=>0.0, :S₂=>0.0] - # ╔═╡ f4d49b1d-9105-4050-a9d4-196fa00a0591 # tspan = missing tspan = (0.0, 150.0) -# ╔═╡ 777ce59f-c849-4a2e-a6dc-ae309d2a2e7c -# parms = missing -parms = [:k=>3.0, :Smax=>150.0, :v=>1e-3, :R=>5.0, :S1res=>10.0, :r=>5] - -# ╔═╡ 43b83336-aea6-4914-bc26-b2e84994ce57 -# oprob = missing -oprob = ODEProblem(sys_irrigation, u0, tspan, parms) - # ╔═╡ 923d04ce-b4d2-44b0-afff-7062c4628ad0 md""" Declare the Turing model. Make sure you take both experiments into account for optimizing $k$ and $S_{max}$. """ # ╔═╡ 481eb8b9-5de2-4f68-b06a-ec18e054c9f5 -# @model function irrigation_inference(t_meas, S1_meas1, S2_meas1, S1_meas2, S2_meas2) +# @model function irrigation_inference(t_meas) # σ_S1 ~ missing # σ_S2 ~ missing # k ~ missing @@ -171,9 +154,9 @@ Declare the Turing model. Make sure you take both experiments into account for o # S1_meas2 ~ missing # S2_meas2 ~ missing # end -@model function irrigation_inference(t_meas, S1_meas1, S2_meas1, S1_meas2, S2_meas2) - σ_S1 ~ InverseGamma() - σ_S2 ~ InverseGamma() +@model function irrigation_inference(t_meas) + σ_S1 ~ Exponential() + σ_S2 ~ Exponential() # k ~ Uniform(0, 10) k ~ LogNormal() Smax ~ Uniform(100, 200) @@ -181,15 +164,15 @@ Declare the Turing model. Make sure you take both experiments into account for o # For experiment 1: u0 = [:S₁=>0.0, :S₂=>0.0] oprob = ODEProblem(sys_irrigation, u0, tspan, parms) - osol1 = solve(oprob, Tsit5(), saveat=t_meas) - S1_meas1 ~ MvNormal(osol1[:S₁], σ_S1^2 * I) - S2_meas1 ~ MvNormal(osol1[:S₂], σ_S2^2 * I) + osol1 = solve(oprob, AutoTsit5(Rosenbrock23()), saveat=t_meas) + S1_ex1 ~ MvNormal(osol1[:S₁], σ_S1) + S2_ex1 ~ MvNormal(osol1[:S₂], σ_S2) # For experiment 2: u0 = [:S₁=>140.0, :S₂=>135.0] oprob = ODEProblem(sys_irrigation, u0, tspan, parms) - osol2 = solve(oprob, Tsit5(), saveat=t_meas) - S1_meas2 ~ MvNormal(osol2[:S₁], σ_S1^2 * I) - S2_meas2 ~ MvNormal(osol2[:S₂], σ_S2^2 * I) + osol2 = solve(oprob, AutoTsit5(Rosenbrock23()), saveat=t_meas) + S1_ex2 ~ MvNormal(osol2[:S₁], σ_S1) + S2_ex2 ~ MvNormal(osol2[:S₂], σ_S2) end # ╔═╡ df933ae8-1f51-4467-93a7-33f153e5e4f8 @@ -199,7 +182,7 @@ Provide the measurements to the Turing model. # ╔═╡ 0e2aa675-9e09-4e06-b5f8-118707ee652a # irrigation_inf = missing -irrigation_inf = irrigation_inference(t_meas, S1_meas1, S2_meas1, S1_meas2, S2_meas2) +irrigation_inf = irrigation_inference(t_meas) | (S1_ex1 = S1_meas1, S2_ex1 = S2_meas1, S1_ex2 = S1_meas2, S2_ex2 = S2_meas2) # ╔═╡ f7f47956-7c3b-44cc-bff7-fb7d32af874a md""" @@ -300,7 +283,7 @@ osol2_opt = solve(oprob2_opt, Tsit5(), saveat=0.5) # end begin plot(osol2_opt, labels=["S1 sim" "S2 sim"], xlabel="t") - scatter!(t_meas, S1_meas2, label="S1 meas2", color=:blue, title="Experment 2") + scatter!(t_meas, S1_meas2, label="S1 meas2", color=:blue, title="Experiment 2") scatter!(t_meas, S2_meas2, label="S2 meas2", color=:red) end @@ -331,11 +314,7 @@ end # ╠═cef4b9a8-b5bf-4a2d-8a8f-5d8f85534859 # ╠═fc2cabd7-e778-4211-bf87-b5c11ca054c9 # ╟─c0b2db7b-0632-4008-9cff-d5fbf3e59807 -# ╟─d35bbe54-5ebc-4ea5-a6c8-a6419476ec4c -# ╠═10057510-e4b6-4a3e-9d3f-f05effc88a58 # ╠═f4d49b1d-9105-4050-a9d4-196fa00a0591 -# ╠═777ce59f-c849-4a2e-a6dc-ae309d2a2e7c -# ╠═43b83336-aea6-4914-bc26-b2e84994ce57 # ╟─923d04ce-b4d2-44b0-afff-7062c4628ad0 # ╠═481eb8b9-5de2-4f68-b06a-ec18e054c9f5 # ╟─df933ae8-1f51-4467-93a7-33f153e5e4f8 diff --git a/exercises/solved_notebooks/P6_calib/optim_wastewater_treatment_sol.jl b/exercises/solved_notebooks/P6_calib/optim_wastewater_treatment_sol.jl index 2b08ec5..810572b 100644 --- a/exercises/solved_notebooks/P6_calib/optim_wastewater_treatment_sol.jl +++ b/exercises/solved_notebooks/P6_calib/optim_wastewater_treatment_sol.jl @@ -202,7 +202,7 @@ Declare the Turing model function. Sample the flow rate $q$ prior from an unifor """ # ╔═╡ b6bac48a-4a3d-47e4-90ea-788ca20dadff -# @model function wastewater_treatment_inference(C_val) +# @model function wastewater_treatment_inference() # q ~ missing # u0 = missing # tspan = missing @@ -211,7 +211,7 @@ Declare the Turing model function. Sample the flow rate $q$ prior from an unifor # osol = missing # C_val ~ missing # end -@model function ww_treat_inf(C_val) +@model function ww_treat_inf() q ~ Uniform(0, 5.0) u0 = [:C=>3.0, :X=>0.5] tspan = (0.0, 72.0) # the time interval to solve on @@ -219,8 +219,8 @@ Declare the Turing model function. Sample the flow rate $q$ prior from an unifor oprob = ODEProblem(sys_ww_treat, u0, tspan, parms) # or, since initial conditions don't change: # oprob = ODEProblem(sys_ww_treat, [], tspan, parms) - osol = solve(oprob, Tsit5(), saveat=0.1) - C_val ~ Normal(osol[:C][end], 1e-3) + osol = solve(oprob, AutoTsit5(Rosenbrock23()), saveat=0.1) + C ~ Normal(osol[:C][end], 1e-3) end # ╔═╡ b3a40556-0c00-4f6d-8cd9-c5fca79d8bbf @@ -232,6 +232,12 @@ Define the desired value for the organic waste with the variable name `C_val`. # missing C_val = 0.28 +# ╔═╡ b828105e-2677-4335-ba93-3bfd2b117220 +md"Instantiate the Turing model and condition it with the observed value of $C$" + +# ╔═╡ 97aca933-ebed-4b97-a872-a3d8ad533130 +ww_model = ww_treat_inf() | (C = C_val,) + # ╔═╡ ee1ffc12-55a1-47ef-ac5b-33148706a09b md""" Provide the measurement to the Turing model and optimize the prior. Do this with `MLE` method and Nelder-Mead. Store the optimization results in `results_mle`. @@ -239,7 +245,7 @@ Provide the measurement to the Turing model and optimize the prior. Do this with # ╔═╡ afc035be-075b-464b-8ba2-20235082f005 # results_mle = missing -results_mle = optimize(ww_treat_inf(C_val), MLE(), NelderMead()) +results_mle = optimize(ww_model, MLE(), NelderMead()) # ╔═╡ 3ee8121e-3e78-4901-a32d-f04d0c6a0996 md""" @@ -338,6 +344,8 @@ Draw your conclusion. # ╠═b6bac48a-4a3d-47e4-90ea-788ca20dadff # ╟─b3a40556-0c00-4f6d-8cd9-c5fca79d8bbf # ╠═2df409ef-bd95-4ac3-a2b8-c5e17c490eba +# ╟─b828105e-2677-4335-ba93-3bfd2b117220 +# ╠═97aca933-ebed-4b97-a872-a3d8ad533130 # ╟─ee1ffc12-55a1-47ef-ac5b-33148706a09b # ╠═afc035be-075b-464b-8ba2-20235082f005 # ╟─3ee8121e-3e78-4901-a32d-f04d0c6a0996 diff --git a/exercises/solved_notebooks/P7_sens_uncert/sens_intro.jl b/exercises/solved_notebooks/P7_sens_uncert/sens_intro.jl index 28c2b68..3fa7fac 100644 --- a/exercises/solved_notebooks/P7_sens_uncert/sens_intro.jl +++ b/exercises/solved_notebooks/P7_sens_uncert/sens_intro.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.19.46 +# v0.20.21 using Markdown using InteractiveUtils diff --git a/exercises/solved_notebooks/P8_modselect/model_selection_friction_sol.jl b/exercises/solved_notebooks/P8_modselect/model_selection_friction_sol.jl index b6c421a..a311ac0 100644 --- a/exercises/solved_notebooks/P8_modselect/model_selection_friction_sol.jl +++ b/exercises/solved_notebooks/P8_modselect/model_selection_friction_sol.jl @@ -216,20 +216,14 @@ We create a Turing model with both parameters $k$ and $g$, where only the first # ╔═╡ fa3383e4-f5fa-4d47-b21b-684c6aacbd46 @model function linear(t) - σ_x ~ InverseGamma() + σ_x ~ Exponential() k ~ LogNormal() params_linear = [:k => k, :g => g] oprob_linear = ODEProblem(model_linear, u0, tspan, params_linear) osol_linear = solve(oprob_linear, AutoTsit5(Rosenbrock23()), saveat=t) - x ~ MvNormal(osol_linear[:x], σ_x^2 * I) + x ~ MvNormal(osol_linear[:x], σ_x) end -# ╔═╡ b6aada44-ef32-4c99-9939-47eedb3fc235 -md""" -!!! note - We use the auto-switching `AutoTsit5(Rosenbrock23())` solver here as the calibration of this model is otherwise quite unstable and will fail occassionally (it's not a problem if you didn't do this). -""" - # ╔═╡ 1f36c90e-ec20-4bdb-8bd0-5f44f03a1748 mod_linear_cond = linear(t) | (x = d,) @@ -457,11 +451,11 @@ Let's now estimate the distance assuming a noisy error and no calibrated paramet # ╔═╡ 965ebc6f-b440-42da-8aa7-00b22e6d0a3f @model function no_resistance(t) - σ_x ~ InverseGamma() + σ_x ~ Exponential() params_nr = [:g => g] oprob_nr = ODEProblem(model_nr, u0, tspan, params_nr) osol_nr = solve(oprob_nr, AutoTsit5(Rosenbrock23()), saveat=t) - x ~ MvNormal(osol_nr[:x], σ_x^2 * I) + x ~ MvNormal(osol_nr[:x], σ_x) end # ╔═╡ d7514616-736b-4cc3-a562-c23fb3d6924d @@ -574,12 +568,12 @@ sol_quad = solve(prob_quad); # ╔═╡ b0908503-1523-4a35-bc2f-5087650e537e @model function quadratic1(t) - σ_x ~ InverseGamma() + σ_x ~ Exponential() k ~ LogNormal() params_quad = [:k => k, :g => g] oprob_quad = ODEProblem(model_quad, u0, tspan, params_quad) osol_quad = solve(oprob_quad, AutoTsit5(Rosenbrock23()), saveat=t) - x ~ MvNormal(osol_quad[:x], σ_x^2 * I) + x ~ MvNormal(osol_quad[:x], σ_x) end # ╔═╡ 8d95ddb9-bd58-4aea-9c8f-28e2f2119be2 @@ -673,13 +667,13 @@ end # ╔═╡ ba0a576e-4b9b-49e1-bfbc-33955dcec808 @model function quadratic2(t) - σ_x ~ InverseGamma() + σ_x ~ Exponential() k₁ ~ LogNormal() k₂ ~ LogNormal() params_quad2 = [:k₁ => k₁, :k₂ => k₂, :g => g] oprob_quad2 = ODEProblem(model_quad2, u0, tspan, params_quad2) osol_quad2 = solve(oprob_quad2, AutoTsit5(Rosenbrock23()), saveat=t) - x ~ MvNormal(osol_quad2[:x], σ_x^2 * I) + x ~ MvNormal(osol_quad2[:x], σ_x) end # ╔═╡ edd1c981-56f6-44d8-a805-0b85c963b39b @@ -772,13 +766,13 @@ sol_power = solve(prob_power); # ╔═╡ 4b53805c-c9ed-4708-b8a1-4922ce6d5a5d @model function power(t) - σ_x ~ InverseGamma() + σ_x ~ Exponential() k ~ LogNormal() r ~ LogNormal() params_power = [:k => k, :r => r, :g => g] oprob_power = ODEProblem(model_power, u0, tspan, params_power) osol_power = solve(oprob_power, AutoTsit5(Rosenbrock23()), saveat=t) - x ~ MvNormal(osol_power[:x], σ_x^2 * I) + x ~ MvNormal(osol_power[:x], σ_x) end # ╔═╡ 4e579e09-2476-432d-8eb4-7a1080376998 @@ -994,7 +988,6 @@ md""" # ╟─85471dd1-df76-4718-a7ee-d501bb8a2030 # ╟─4bde3bb3-5dcd-42d4-ae48-b922e99a04bc # ╠═fa3383e4-f5fa-4d47-b21b-684c6aacbd46 -# ╟─b6aada44-ef32-4c99-9939-47eedb3fc235 # ╠═1f36c90e-ec20-4bdb-8bd0-5f44f03a1748 # ╠═c4de5cb6-4b84-46b0-acc5-01aea794de3e # ╠═3153b018-57f8-4d7a-bfb9-c4ca2886363a diff --git a/exercises/solved_notebooks/P8_modselect/model_selection_intro.jl b/exercises/solved_notebooks/P8_modselect/model_selection_intro.jl index 5b42e6a..d7471d0 100644 --- a/exercises/solved_notebooks/P8_modselect/model_selection_intro.jl +++ b/exercises/solved_notebooks/P8_modselect/model_selection_intro.jl @@ -173,25 +173,14 @@ osys_log = convert(ODESystem, growth_log) # ╔═╡ cd8a7ba1-194a-4b42-8c87-2c9b1fe6b475 md""" -Setting initial conditions, timespan and parameter values: +Setting the initial conditions and parameter values. """ # ╔═╡ cbb2ca49-b019-495b-9310-83fcc00cad26 u0_log = [:W => 2.0] -# ╔═╡ be565a3c-31b6-4df1-b73b-08f308a8c09b -md""" -For the sake of clarity, we will use the variables `μ_log` and `Wf_log` to store the parameter values. -""" - -# ╔═╡ f62806d1-77e1-470b-9711-33a924c788cc -μ_log = 0.07 - -# ╔═╡ 546ed163-26a7-4235-982f-7568ed609488 -Wf_log = 10.0 - # ╔═╡ 0da53fa2-5a42-46e6-8bd8-45d6aa903d46 -params_log = [:μ => μ_log, :Wf => Wf_log] +params_log = [:μ => 0.07, :Wf => 10.0] # ╔═╡ b1298f40-4696-49d0-ac94-896e0cdbc996 md""" @@ -200,7 +189,7 @@ Creating and solving the ODEProblem and plotting results: # ╔═╡ 270647d2-1371-4272-8bc1-3a6ad77bc716 oprob_log = ODEProblem(growth_log, u0_log, tspan, params_log); -# Also possible here if initial conditions and parameter values are defined in the catalyst model: +# Also possible here since the initial conditions and parameter values are defined in the catalyst model: # oprob_log = ODEProblem(growth_mod_log, [], tspan, []) # ╔═╡ ac235d86-1d93-4944-aa89-1b4fd38f0e6e @@ -215,7 +204,7 @@ end # ╔═╡ e2955af0-edf0-4702-9a8f-478141ffdc3b md""" -We can see that the model does not predict well the data set for the considered parameter values. Thus we will use the data to both calibrate the model parameters and assess the quality of the fit. +We can see that the model does not predict the data well using these parameter values. Thus we will use the data to both calibrate the model parameters and assess the quality of the fit. """ # ╔═╡ c73669c4-d7af-4877-b32d-6499f339e27a @@ -230,22 +219,21 @@ We declare our Turing model function: # ╔═╡ 169a67ff-55fb-4d1d-b98b-126f4af47e77 @model function growth_log_fun(t_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - Wf ~ LogNormal() + σ_W ~ Exponential() + W0 ~ LogNormal(log(1)) + μ ~ LogNormal(log(0.1)) + Wf ~ truncated(LogNormal(log(10)), lower = 0.1) # prevent zeros in denominator u0_log = [:W => W0] - params_log = [:μ => μ, :Wf => Wf] - oprob_log = ODEProblem(growth_log, u0_log, tspan, params_log) + parms_log = [:μ => μ, :Wf => Wf] + oprob_log = ODEProblem(growth_log, u0_log, tspan, parms_log) osol_log = solve(oprob_log, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_s ~ MvNormal(osol_log[:W], σ_W^2 * I) - return osol_log # optionally, to be used with MCMC + W_s ~ MvNormal(osol_log[:W], σ_W) end # ╔═╡ bfb26b9d-6969-457c-8567-03422a7f2a93 md""" -!!! note - We use the auto-switching `AutoTsit5(Rosenbrock23())` solver here as the calibration of this model is otherwise quite unstable and will fail occassionally (it's not a problem if you didn't do this). +!!! tip + This model can change between being non-stiff and stiff based on the sampled parameter values. You can use an auto-switching solver such as `AutoTsit5(Rosenbrock23())` here to make calibration more stable. """ # ╔═╡ aa3e553d-2731-42d6-b0e6-1821e4d7f4d4 @@ -574,15 +562,15 @@ Declare the Turing model function. # ╔═╡ 89bf91c6-117c-4647-bfb9-6fc9b8dcfb5f @model function growth_exp_fun(t_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - Wf ~ LogNormal() + σ_W ~ Exponential() + W0 ~ LogNormal(log(1)) + μ ~ LogNormal(log(0.1)) + Wf ~ LogNormal(log(10)) u0_exp = [:W => W0] - params_exp = [:μ => μ, :Wf => Wf] - oprob_exp = ODEProblem(growth_exp, u0_exp, tspan, params_exp) + parms_exp = [:μ => μ, :Wf => Wf] + oprob_exp = ODEProblem(growth_exp, u0_exp, tspan, parms_exp) osol_exp = solve(oprob_exp, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_s ~ MvNormal(osol_exp[:W], σ_W^2 * I) + W_s ~ MvNormal(osol_exp[:W], σ_W) end # ╔═╡ b9a8c0cb-fe38-448e-aa4c-e9422f24a4a4 @@ -746,7 +734,7 @@ $W_0$ = 2.0, $\mu$ = 0.09 and $D$ = 0.04. # ╔═╡ 68895842-0a1b-4236-b159-19a2805ea44d growth_gom = @reaction_network begin - μ-D*log(W), W --> 2W + μ-d*log(W), W --> 2W end # ╔═╡ ec12d3c4-6b48-4f3a-b463-0fc67d3adff3 @@ -757,15 +745,15 @@ Declare the Turing model. Take the same priors as before. # ╔═╡ 1f8b10ba-f380-492e-800c-2673774cb25c @model function growth_gom_fun(t_meas) - σ_W ~ InverseGamma() - W0 ~ LogNormal() - μ ~ LogNormal() - D ~ LogNormal() + σ_W ~ Exponential() + W0 ~ truncated(LogNormal(log(1.0)), lower = 1e-5) # prevent log(0) + μ ~ truncated(LogNormal(log(0.1)), upper = 1.0) # prevent overly large exponential growth from crashing the solver - this model has no carrying capacity to slow things down + d ~ LogNormal(log(0.1)) u0_gom = [:W => W0] - params_gom = [:μ => μ, :D => D] - oprob_gom = ODEProblem(growth_gom, u0_gom, tspan, params_gom) + parms_gom = [:μ => μ, :d => d] + oprob_gom = ODEProblem(growth_gom, u0_gom, tspan, parms_gom) osol_gom = solve(oprob_gom, AutoTsit5(Rosenbrock23()), saveat=t_meas) - W_s ~ MvNormal(osol_gom[:W], σ_W^2 * I) + W_s ~ MvNormal(osol_gom[:W], σ_W) end # ╔═╡ 82bee644-7b9a-42fa-a32f-76b3eb5038b5 @@ -778,11 +766,11 @@ growth_gom_cond_mod = growth_gom_fun(t_meas) | (W_s = W_meas,) # ╔═╡ a438470b-e037-4b65-b6b0-fb9c202e74a4 md""" -Optimize the priors ($\sigma_W$, $W_0$, $\mu$ and $D$). Do this now with `MAP` method and Nelder-Mead. Store the optimization results in `results_gom_map`. +Optimize the priors ($\sigma_W$, $W_0$, $\mu$ and $D$). Do this now with the `MAP` method and Nelder-Mead. Store the optimization results in `results_gom_map`. """ # ╔═╡ f4c3dcd0-70ec-49ab-bbad-04b550da481c -results_gom_mle = optimize(growth_gom_cond_mod, MLE(), NelderMead()) +results_gom_map = optimize(growth_gom_cond_mod, MAP(), NelderMead()) # ╔═╡ df166b45-3c28-4bad-95b8-f43429a58fd0 md""" @@ -790,7 +778,7 @@ Visualize a summary of the optimized parameters. """ # ╔═╡ a2738236-86a1-419a-9617-b1229f7c9240 -coeftable(results_gom_mle) +# coeftable(results_gom_map) # ╔═╡ 59f78a22-cac1-49a7-b3e6-9d74b694be64 md""" @@ -798,13 +786,13 @@ Get the optimized values and assign them to `W0_opt_gom`, `μ_opt_gom` and `D_op """ # ╔═╡ 7c351d72-c5ce-4e5a-b2e5-87ff9dbc70e5 -W0_opt_gom = coef(results_gom_mle)[:W0] +W0_opt_gom = coef(results_gom_map)[:W0] # ╔═╡ d7312308-7c95-4c14-9bd6-56a7ac99d49a -μ_opt_gom = coef(results_gom_mle)[:μ] +μ_opt_gom = coef(results_gom_map)[:μ] # ╔═╡ 8f8b23b3-2e33-46e4-917a-294934fa090e -D_opt_gom = coef(results_gom_mle)[:D] +d_opt_gom = coef(results_gom_map)[:d] # ╔═╡ 4ff1079f-8dfb-4069-9a98-0d263eba9920 md""" @@ -825,7 +813,7 @@ Set up parameter values with optimized parameter values: """ # ╔═╡ 0d225f18-bdf6-49d1-a827-2bc9551c158d -params_opt_gom = [:μ => μ_opt_gom, :D => D_opt_gom] +params_opt_gom = [:μ => μ_opt_gom, :d => d_opt_gom] # ╔═╡ 5d27e199-ad09-454a-8966-1531371ce67b md""" @@ -1014,7 +1002,7 @@ We can repeat the calibration and take into account the priors to obtain the MAP results_log_map = optimize(growth_log_cond_mod, MAP(), NelderMead()) # ╔═╡ 3d5e9bb7-5e65-4de5-a521-6b1f1cb872b3 -coeftable(results_log_map) +# coeftable(results_log_map) # ╔═╡ 3a3cd103-78b1-481a-8746-749d402346d1 md""" @@ -1095,6 +1083,9 @@ WAIC_exp = -2*(lppd_exp - pWAIC_exp) # ╔═╡ 8330fce5-13c9-447e-8368-01b7c7ebc68f results_gom_nuts = sample(growth_gom_cond_mod, NUTS(), N) +# ╔═╡ a907bd5d-455d-4471-94c7-a6b26770f81f +plot(results_gom_nuts) + # ╔═╡ 0d2cab0c-04a2-4c8b-bb4f-2115ae0cf871 lppd_gom = sum(results_gom_nuts.value[:, :lp]) @@ -1158,9 +1149,6 @@ md""" # ╠═ba181db8-d176-4b83-9168-d5939ffe9661 # ╟─cd8a7ba1-194a-4b42-8c87-2c9b1fe6b475 # ╠═cbb2ca49-b019-495b-9310-83fcc00cad26 -# ╟─be565a3c-31b6-4df1-b73b-08f308a8c09b -# ╠═f62806d1-77e1-470b-9711-33a924c788cc -# ╠═546ed163-26a7-4235-982f-7568ed609488 # ╠═0da53fa2-5a42-46e6-8bd8-45d6aa903d46 # ╟─b1298f40-4696-49d0-ac94-896e0cdbc996 # ╠═270647d2-1371-4272-8bc1-3a6ad77bc716 @@ -1347,6 +1335,7 @@ md""" # ╠═eb9fe668-0a0f-4a8f-a793-ba545cb02e85 # ╠═ddf9ae6d-0abe-4c8c-9646-024dd9a41e02 # ╠═8330fce5-13c9-447e-8368-01b7c7ebc68f +# ╠═a907bd5d-455d-4471-94c7-a6b26770f81f # ╠═0d2cab0c-04a2-4c8b-bb4f-2115ae0cf871 # ╠═6900e703-de7d-4ea3-8ac5-0b38b21b083d # ╠═b353eab2-d78f-4c5e-8e38-bd57813b6816 diff --git a/exercises/student_notebooks/P6_calib/calib_irrigation.jl b/exercises/student_notebooks/P6_calib/calib_irrigation.jl index 94fb8d1..f40af3d 100644 --- a/exercises/student_notebooks/P6_calib/calib_irrigation.jl +++ b/exercises/student_notebooks/P6_calib/calib_irrigation.jl @@ -1,5 +1,5 @@ ### A Pluto.jl notebook ### -# v0.20.4 +# v0.20.21 using Markdown using InteractiveUtils